Estimation of activity coefficients for aqueous organic redox flow batteries: Theoretical basis and equations.

The field of aqueous organic redox flow batteries (AORFBs) has been developing fast in recent years, and many chemistries are starting to emerge as serious contenders for grid-scale storage. The industrial development of these systems would greatly benefit from accurate physics-based models, allowing to optimize battery operation and design. Many authors in the field of flow battery modeling have brought evidence that the dilute solution hypothesis (the assumption that aqueous electrolytes behave ideally) does not hold for these systems and that calculating cell voltage or chemical potentials through concentrations rather than activities, while serviceable, may become insufficient when greater accuracy is required. This article aims to provide the theoretical basis for calculating activity coefficients of aqueous organic electrolytes used in AORFBs to provide tools to predict the concentrated behavior of aqueous electrolytes, thereby improving the accuracy of physics-based models for flow batteries.

Introduction

Among energy storage technologies, batteries have recently seen a surge in popularity, owing to their flexibility of installation and lower capital costs required compared with larger centralized storage technologies, allowing the transition from a centralized grid model to that of a network of smart, decentralized microgrids (Gür, 2018). In this context, aqueous organic redox flow batteries (AORFBs) present a number of significant advantages, including: absence of metal-based active materials, leading to a greener and more resilient supply chain of battery materials when compared to metal-based technologies; decorrelation between energy and power sizing, leading to a greater flexibility in design to meet local grid requirements; and finally non-flammability because of the large volumes of water-based electrolytes, leading to a greater reliability compared to technologies presenting the risk of thermal runaway (lithium-ion) or high flammability (hydrogen).

Flow batteries and concentrated solution theory

Since the first modeling studies conducted by NASA (Fedkiw and Watts, 1984), computational models have become increasingly important for optimizing operating and design parameters of the batteries. Although it is common practice to assume the aqueous electrolyte solutions flowing in the battery to be ideal, i.e., where the ions are assumed not to interact with each other (the dilute solution hypothesis), a growing number of authors are opposing this hypothesis, advocating for using concentrated solution theory instead, which accounts for non-idealities of mixing. In order to assess the current state of the literature consensus on the topic, we conducted a meta-analysis on a literature sample (N = 66) focusing on the development of computational models for redox-flow batteries. Based on the featured equations and their hypotheses and justifications, papers were classified in four categories.

In category A, equations featured uncorrected concentrations (e.g., in voltage predictions), but the dilute solution hypothesis was not mentioned explicitly. As an illustration, You et al. (2009) manipulated concentrations to calculate the open-circuit potential (OCP) of an electrode aswhere is the reference thermodynamic potential (V vs. SHE), the ideal gas constant, the absolute temperature, Faraday’s constant and and (mol/L) the concentrations of reduced and oxidized species, respectively, of the redox pair in the considered electrode. You et al. (2009) noted that “a total value of 140 mV was added to the simulated results in order to match the experimental data” but the fact that this voltage discrepancy may arise from the hypothesis of ideal solutions was not mentioned.

In category B, the dilute solution hypothesis was mentioned explicitly: For example, Shah et al. (2011) mention that “concentrated-solution theory is straightforward but […] several parameters are not known for the present system” so that as in You et al., “A total value of 131 mV was [added] from fitting the experimental data to simulation results”.

In category C, authors either attempted to justify the dilute solution hypothesis, or at least question its validity. Knehr and Kumbur (2011) belong to the former category, with very debatable justifications e.g., “in a VRFB, the liquid electrolytes are typically well-circulated, yielding negligible interactions among the ions”. They did rightly, however, introduce an additional term in the open-circuit potential to account for possible imbalances in the concentration of supporting ionswhere corresponds to the concentration of species and to their stoichiometric coefficient in the half-cell reaction, including exchanged ions (protons in the case of Knehr and Kumbur (2011)). Stephenson et al. (2012) instead acknowledged that standard potentials in the Nernst equation may be shifted from non-ideal effects “Like Shah et al.[…] we also observed a shift in equilibrium potential. It may be attributed to the non-ideal solution, in which the activity coefficient is not equal or close to one.” leading to the first mention of the possibility that activity coefficients may deviate from their ideal value in flow batteries since Heintz and Illenberger (1998). No attempt was made, however, at quantifying these activity coefficients and the discrepancies arising from the behavior of non-ideal solutions.

In category D, authors attempted to manipulate equations in accordance with concentrated solution theory and sometimes advocated against the dilute solution hypothesis. Hayer and Kohns (2020) for example state that “fitting kinetic parameters to measurement data of voltage losses can lead to ambiguous results if only the idealized OCV, obtained by neglecting the activity coefficients, is used in that evaluation” and therefore re-introduce the more fundamental form of the Nernst equation for the calculation of the open-circuit voltagewhere corresponds the chemical activity of species in the electrodes, a quantity central to concentrated solution theory, which will be the focus of the present publication.

It seems apparent from the visualization of the dilute versus concentrated debate on Figure 1 that although the dilute solution hypothesis was very prevalent between 2000 and 2010 (category A + B), a few papers started either questioning the assumptions behind it, or at least acknowledging that the dilute solution hypothesis should be justified (category C), between 2010 and 2018. Although a few papers re-introduced the notion of chemical activity (which was present in the paper from Heintz and Illenberger (1998) but considered negligible in the following decade), a vast majority of the literature preferred the use of concentration. Between 2018 and today, however, it seems that the literature reviewed in this analysis has polarized into two main schools of thoughts: On the one hand, papers using concentrations in their parametrization, not mentioning the dilute solution hypothesis and its implications, and on the other hand authors advocating for the manipulation of equations in line with concentrated solution theory.

Figure 1

Visual representation of the dilute vs. concentrated debate in flow battery modeling literature

Papers are ordered in 4 categories. (A) Papers using the dilute solution hypothesis without mentioning it; 2000–2010: [6] Sukkar and Skyllas-Kazacos (2003), [7] Li and Hikihara (2008), [8] You et al. (2009); 2010–2015: [11] Tang et al. (2011), [12] Qiu et al. (2012), [13] Tang et al. (2012), [14] Badrinarayanan et al. (2014); 2015-2018: [26] Zhou et al. (2016), [27] Nikonenko et al. (2016), [28] Yan et al. (2016), [29] Milshtein et al. (2017), [30] Bhattacharjee and Saha (2017); 2018-2022: [45] Lei et al. (2018), [46] Murali et al. (2018), [47] Chen et al. (2019) [48] Xu et al. (2020) [49] Chen et al. (2021a), [50] Chen et al. (2021b) [51] Howard Amanda and Tartakovsky (2021).

(B) Papers mentioning the dilute solution hypothesis but not justifying it; 2000–2010: [9] Shah et al. (2008), [10] Al-Fetlawi et al. (2009); 2010–2015: [15] Shah et al. (2011), [16] Vynnycky (2011), [17] Knehr et al. (2012), [18] Yin et al. (2014) [19] Chen et al. (2014); 2015-2018: [31] Xu and Zhao (2015), [32] Yang et al. (2015), [33] Boettcher et al. (2015), [34] Zhou et al. (2015), [35] Lei (2015), [36] Gandomi et al. (2016), [37] Chu et al. (2016), [38] Chen et al. (2017), [39] Elgammal et al. (2017); 2018-2022: [52] Ronen et al. (2018), [53] Oh et al. (2019), [54] Delgado et al. (2020), [55] Cacciuttolo et al. (2021).

(C) Papers justifying or discussing the dilute solution hypothesis; 2010–2015: [20] Knehr and Kumbur (2011), [21] Stephenson et al. (2012), [22] Wandschneider et al. (2014), [23] Hudak (2014), [24] Wei et al. (2014); 2015-2018: [40] Darling et al. (2016), [41] Darling et al. (2016), [42] Moshtarikhah et al. (2017); 2018-2022 [56] Murthy et al. (2018) [57] Al-Yasiri (2020).

(D) Papers using equations coherent with concentrated solutions theory; 2010-2015: [25] Ashraf Gandomi et al. (2014); 2015–2018: [43] Pavelka et al. (2015), [44] Xu et al. (2017); 2018-2022: [58] Intan et al. (2018), [59] Sijabat et al. (2019), [60] Hao et al. (2019), [61] Kleinsteinberg et al. (2020), [62] Hayer and Kohns (2020), [63] del Olmo et al. (2020), [64] Crothers et al. (2020) [65] Vardner et al. (2020) [66] Mourouga et al. (2022). Distance from the center and clockwise angle are proportional the publication date. The 5 earliest publications (Reid and Gahn (1977), Fedkiw and Watts (1984), Zawodzinski et al. (1993), Mohammadi et al. (1997), Heintz and Illenberger (1998)) are shown for reference but were not categorised.

Activity coefficients impact the OCV of the cell and typically introduce a few tens of mV of discrepancy with the same equations parametrized with concentrations (Hayer and Kohns, 2020; del Olmo et al., 2021; Mourouga et al., 2022). Although this value may sound negligible for performance models over a single cycle, modeling a flow battery over several thousands of charge-discharge cycles requires accurate transport models. Transport models, in turn, are highly impacted by the chemical potentials of species mixing in the battery, which are directly impacted by thermodynamic activity according to the formulawhere is the chemical potential, is a reference chemical potential in standard conditions and is the thermodynamic activity, capturing interactions in highly concentrated solutions.

In addition to chemical potentials, transport models are also highly impacted by cycling times which, in turn, correspond to the time at which the cell voltage reaches its cut-off values. Therefore, uncertainties of a few tens of mV on the cell voltage will result in a few minutes of uncertainty on the first cycling time, which will quickly propagate to several hours over thousands of cycles. Accurate OCV calculations, using activity estimation as in Equation (3) rather than concentrations as in Equation (2), thus significantly improve the accuracy of flow battery models over long time scales.

Furthermore, thermodynamic activity also impacts other characteristics of the electrolytes, such as their conductivity and heat capacities. Accurate physics-based shunt currents and thermal models therefore also greatly benefit from accurate activity coefficient values over a wide range of composition and temperature.

To the best of our knowledge, no publication has so far addressed the problem of measuring activity coefficients of flow battery electrolytes, especially for aqueous organic redox-flow batteries. The estimation of activity coefficients for flow batteries will therefore be the focus of this series of publication, in which the present one will introduce the theoretical foundation for the calculation of activity coefficients from experimental data on the properties of solutions.

A short history of activity

The term chemical activity was first coined by Lewis (1907) at a time where two different approaches were competing for predicting observable properties of aqueous electrolytes such as vapor pressure, osmotic pressure, and freezing point. The first one was led by Van’t Hoff, Nernst, and Arrhenius and aimed to apply laws demonstrated from the kinetic theory of ideal gases to chemical systems, whereas the second one, led by Gibbs (1879), Duhem, Planck and Lewis aimed to apply the fundamental laws of physics and thermodynamics to chemical systems.

The main difference between the two approaches lies in how they express the chemical potential of solutes dissolved in a solvent: The kinetic theory of ideal solutions approach relies on the assumption that charged particles in an aqueous solution behave similarly to ideal gases, which allows to express the chemical potential as a function of the logarithm of the molality (in moles of solute per kilogram of solvent)where is a chemical potential taken in standard conditions of pressure and temperature, and is a standard molality, chosen as 1 mol kg−1 for convenience.

The approach championed by Lewis, on the other hand, defines the notions of chemical activity , which is equal to the molality multiplied by an activity coefficient accounting for interactions between the mixing species

It should be noted that the standard state and standard molality have the same physical meaning in both frameworks: they represent the hypothetical case of an ideal (non-interacting) mixture at a molality of one mole of solute per kilogram of solvent. Any calculated chemical potential is therefore a deviation from this ideal state, which follows from the fact that only differences in Gibbs free energy between two thermodynamic states, as opposed to absolute Gibbs free energies , are measurable.

Debye and Hückel (1923) provided the first theoretical calculation of activity coefficients, coherent with experimental measurements in dilute solutions. Their main idea was that charged ions obeyed the law of electrostatic interactions, and that through the application of the Boltzmann principle and the Poisson equation, one could derive the changes in potential energy of ionic solutions and thus, predict measurable quantities such as freezing point depression: For all electrolytes, in the limit for low concentrations, the percentage deviation of the freezing point depression from the classical value is proportional to the square root of the concentration. This theorem (and the two that follow, see Debye et al. (1923) for the English translation) constitute the Debye–Hückel limiting law and state that in the limit of low concentrations, activity coefficients can be interpreted as the influence of the electrostatic interaction on electrolyte solutions (Hückel and Braus, 2020).

Scatchard (1936) corrected some assumptions from the Debye–Hückel theory and replaced the use of concentrations (moles of solute per volume of solution) with molalities (moles of solute per mass of solvent), then Guggenheim (1935) and Guggenheim and Prue (1954) demonstrated the relationship between chemical activity and other properties of aqueous electrolytes, including enthalpies of mixing and heat capacity. Guggenheim and Turgeon (1955) further cemented the idea that higher order interactions were in play in more concentrated solutions.

The most general theory was demonstrated by Pitzer (1973), who showed that the virial theorem could be applied to aqueous electrolyte solution, leading to an expansion of the osmotic pressure and the excess Gibbs free energy (and thus, interactions) in powers of the molality . Within this framework, activity coefficients can be viewed as an expression of different interactions in solution: represents overall repulsive (long-range) interactions, and represents overall attractive (higher-order, short-range) interactions. The higher-order interactions become dominant when the average distance between molecules is reduced, which therefore come to dominate in concentrated solutions.

The very extensive work of Pitzer and his successors (Pitzer and Mayorga (1973); Bradley and Pitzer (1979); Clarke and Glew (1966); Clarke and Glew (1980); Clarke and Glew (1985); Craft and Van Hook 1975a; Craft and Van Hook 1975b; Tanner and Lamb (1978); Fortier et al. (1973)) solidified his thermodynamic derivations and validated them with vapor pressure, liquid junction potential and calorimetry measurements on many different aqueous solutions and mixtures. Although this approach allows to derive relationships between properties of a mixture that could appear disconnected (osmotic pressure and excess heat capacity for example), it should be emphasized that it remained a semi-empirical approach, where many phenomenological coefficients had to be fitted to experimental results in order to achieve a satisfying predicting power.

Other mathematical frameworks relying on the calculation of the excess Gibbs free energy of a mixture, but with an emphasis on computation rather than experimentation, emerged after the theory of Pitzer, such as the universal quasi-chemical equation (UNIQUAC) developed by Abrams and Prausnitz (1975) and extended by other authors to UNIFAC (Fredenslund et al., 1975; Maurer and Prausnitz, 1978) or GEQUAC (Egner et al., 1997).

According to Abrams and Prausnitz themselves ”the UNIQUAC equation gives good representation of both vapor-liquid and liquid-liquid equilibria for binary and multicomponent mixtures containing a variety of non-electrolyte components such as hydrocarbons, ketones, esters, amines, alcohols, nitriles, etc …”.

Other computational methods, such as COSMO (Klamt and Schüürmann, 1993) which later branched into COSMO-RS (Klamt, 1995) and COSMOSPACE (Klamt et al., 2002), were developed to bridge macroscopic statistical thermodynamics and quantum chemistry. Quantum chemistry is a fast-developing and exciting field, and future models may hold very real interests for the development of flow batteries, such as screening of new molecules or non-aqueous solvents.

In the present work, however, it was chosen to favor the semi-empirical approach of Pitzer-Clarke to estimate activity coefficients, as the focus of our work was to compute the properties of available flow battery electrolytes from a set of experiments. In the more general case where novel organic electrolytes would be investigated in silico, this approach might prove insufficient to screen the most promising electrolytes among a set of molecules. It is, however, undeniably the most accurate, provided experiments can be readily performed on the aqueous solutions and mixtures (Klamt et al., 2002).

Calculation of activity coefficients

We will use the approach developed by Pitzer (1973) and Clarke and Glew (1966) to demonstrate a combined experimental and computational method for the estimation of activity coefficients from infinitely dilute solution to saturation, and a range of temperature ranging from 0 to 60°C, at atmospheric pressure, from selected experimental data of measurable properties of aqueous solutions.

Ideal solution model

Let us consider an aqueous solution system comprised moles of an XY solute dissolved in moles of water. The solute molality , in moles of solute per kilogram of solvent, is defined aswhere is the molar mass of water expressed in g/mol. Because it is assumed that the XY solute dissolves into cations and anions the individual ionic molalities are defined asand the ionic strength of the system is defined aswhere is the valence of ionic species dissolved in solution and their molality. In the case of an ideal solution, it is assumed that the solute is fully dissociated and that the ions do not interact with each other (through electrostatic or Van Der Waals interactions, for example).

Within this framework, the kinetic theory of ideal gases developed by Maxwell and Boltzmann in the 1870s should be applicable to aqueous solutions, and their colligative properties (the osmotic pressure in Pa, the ionic conductivity in S/m, etc …) should be correlated with the molality .

Experimental deviations from the ideal solution model

The osmotic pressure of aqueous solutions was shown by Van’t Hoff in 1887 to follow the relationshipWhere is the ideal gas constant, the temperature, the molality of solute and the Van’t Hoff number.

Van’t Hoff showed that within an ideal solution framework, should be equal to the number of ions the solute dissociates into , but his experimental measurements pointed towards non-constant values for , below that of the expected value at low concentrations and greater at high concentrations. Similar deviations from ideality would be observed by Arrhenius in 1884 (Hingerl, 2012) on the electrical conductivity of solutions, by Jones and Getman (1903) on their freezing point depression and by Marshall (1906) on their vapor pressure at saturation. Lewis and Randall (1921) challenged the values of the equilibrium constants measured by Guldberg and Waage in their law of mass action reported in 1899: “to those who have not examined this question closely, it may be surprising to learn how far [the equilibrium constant] is from being a constant”. He would be the first to formally introduce the notion of the activity coefficient , a mathematical function lower than 1 at low concentrations and greater than 1 at higher concentrations, to correct the ”apparent molality” of solute and capture the non-linearities observed experimentally on the colligative properties of solution.

Using modern notations, the introduction of the osmotic coefficient (proportional to the activity of water and related to through the Gibbs-Duhem equation) in the law demonstrated by Van’t Hoff, allows to calculate the observed deviation from ideality of the osmotic pressure with a simple expression

In the general case, the Van’t Hoff number is therefore equal to , and in the particular case of ideal solutions where , we obtain the real number of ions the solute dissociates into , as demonstrated by Van’t Hoff in 1887.

Activity, Gibbs free energy of solutions and chemical potential

In 1873, Gibbs introduced the notion of Gibbs free energy of a solution to quantify the amount of useful work which could be extracted from it, and the notion of chemical potential (molar free energy) to quantify how the available energy would evolve as a function of the number of moles of each component. For the aqueous solution previously defined (as demonstrated in A)where is the total number of mole of ions in solution, the chemical potential of the solute, the number of moles of water and its chemical potential.

Lewis would build on the work of Gibbs in 1907 by decomposing the Gibbs free energy of solution as the sum of two terms: the ideal Gibbs free energy , capturing the ideal solution behavior as predicted by the kinetic theory of gases, and the apparent molal excess Gibbs free energy , capturing the observed deviations from ideal solution behavior:where can be demonstrated to be proportional to the activity coefficient and osmotic coefficient according to the relationship (Clarke and Glew, 1985)

More details on the derivation of Equation (12) and Equation (14) can be found in Appendix.

Continuum model of non-ideal solutions

It was shown by Pitzer (1973) that the behavior of ions dissolved in a solvent could be approached by considering interacting charged hard spheres (the ions) in a continuous medium (the solvent). In this framework, the average over time of the total kinetic energy of the ions, bound by long-range electrostatic interactions and short-range van der Waals interactions, was related to that of the total potential energy of the solution (its Gibbs free energy) through the virial theorem. This resulted in a development of the osmotic pressure and the Gibbs free energy of solution in powers of the molality , which would allow to calculate the coefficients and accurately as a function of composition.

Virial expansion of the osmotic pressure

The virial expansion applied to the osmotic coefficient would yield the following dependence (Pitzer and Mayorga, 1973):Where coefficient A is derived from the Debye–Hückel theory and only depends on the solvent, and the following terms (coefficients Q, B, C, D, E, …) are virial coefficients corresponding to solute interactions (e.g., electrostatic or van der Waals), developed in powers of the molality . The Pitzer theory does not provide a framework for the calculation of the virial coefficients, which need to be fitted to experimental data.

From the Gibbs–Duhem equation explicited in Appendix, it can be demonstrated that the dependency of the activity coefficient of solute with the molality follows the relationshipby inserting Equation (15) and Equation (16) in Equation (14), we get the following expression for the apparent molal excess Gibbs free energy

In short, once the Pitzer coefficients B, Q, C, D, … are known, they allow to calculate the variations of the excess Gibbs free energy, the activity coefficient, and the osmotic coefficient as a function of the molality . This, in turn, allow to accurately capture experimental non-linearities observed in the osmotic pressure, the ionic conductivity, the freezing point depression or the vapor pressure of aqueous electrolyte solutions as a function of composition.

However, the values of the coefficients A, B, Q, C, D, …in (Equation 15), (Equation 16), (Equation 17) as reported in Pitzer and Mayorga (1973), are given at 25°C. Although colligative properties of aqueous solutions, such as the ionic conductivity or the osmotic pressure, are not constant with temperature, the Pitzer coefficients should also exhibit variations with the temperature. This dependence was further investigated by a few authors (Cepeda, 1977; Phutela and Pitzer, 1986) and can be linked to the excess enthalpy of solution via a Taylor expansion in the neighborhood of a reference temperature.

Taylor expansion of the excess enthalpy

When mixing two non-ideal solutions of different compositions, an amount of heat, proportional to the enthalpy of mixing, is generated. Similarly, the heat capacity of the solution varies as a function of the composition and the temperature, in a way which is related to the enthalpy of mixing. Both phenomena are fundamentally related to the interactions between ions in solution and therefore, to the excess Gibbs free energy of the mixture.

Clarke and Glew (1966) demonstrated that the enthalpy of mixing (J/mol) and the excess heat capacity (J/mol/K) being well-behaved functions of the temperature, a Taylor expansion in powers of , where is a reference temperature (usually chosen as 25°C), could be applied as:

The Gibbs–Helmholtz equation gives the variations of the activity coefficient as a function of temperature from the enthalpy of mixing:

Therefore, the temperature dependence of the activity and osmotic coefficients is not an expansion in , as for the enthalpy, but rather in .

Furthermore, the relationship demonstrated in Clarke and Glew (1966), explicited in a more complete form in Clarke and Glew (1985) showed that one could introduce the variable to express the integration as:an expression which is explicited in the temperature matrices given in Appendix.

Matrix expression

Expressions for the activity coefficient, the molal enthalpy of mixing and the excess heat capacity of solution as a function of both composition and temperature are better represented in matrix form aswhere the matrices are arrays of molality terms which represent the interactions between the species in solution, explicited through the virial expansion as introduced in (Equation 15), (Equation 16), (Equation 17), and the matrixes are temperature arrays, explicited through the Taylor expansion as introduced in (Equation 18), (Equation 19), (Equation 20). A complete expression of these matrixes can be found in B.

V is a matrix of virial coefficients of the formwhere the coefficients are given by the Debye-Hückel theory as computed by Clarke and Glew (1980) and the coefficients are semi-empirical Pitzer coefficients fitted to available experimental data. This data may include colligative properties of non-ideal aqueous solutions such as vapor pressure, enthalpy of mixing, heat capacity or freezing point depression, as a function of composition and temperature. The fitting procedure and its associated experimental cost will be the focus of the next sections.

Application to reference aqueous electrolytes

Complete virial matrix for NaCl

Clarke and Glew (1985) compiled experimental data of various thermodynamic properties (activity coefficient, osmotic coefficient, solubility, heat of dilution, heat of solution, heat capacity) for aqueous sodium chloride solutions at different composition and temperatures from more than 100 references. Details on the treatment of experimental data from different techniques (isopiestic measurements, freezing point depression measurement, electromotive force of concentration cells, flow calorimetry) and their respective weighing are reported in Clarke and Glew (1985).

The result is the following matrix of virial coefficientswhere the coefficients are given by the Debye–Hückel theory as computed by Clarke and Glew (1980) and the coefficients are Pitzer coefficients which values can be found in Clarke and Glew (1985). In the present work these values, as well as the theoretical expressions from Eq.(B.1)-(B.4) (in Appendix) were entered in a MATLAB R2020b function in order to evaluate the osmotic and activity coefficients, the excess enthalpy and heat capacity as a function of composition and temperature.

The validation dataset, taken from reference literature data, was compiled from Messikomer and Wood (1975) for the excess enthalpy, from Tanner and Lamb (1978) for the excess heat capacity and from Gibbard et al. (1974) for the osmotic coefficient. The resulting surface plots are given on Figure 2, which shows that the complete Virial matrix calculated by Clarke and Glew (1985) allows to accurately estimate properties of aqueous sodium chloride as reported by other authors.

Figure 2

Graphs showing the agreement between (A) excess enthalpy (B) excess heat capacity (C) osmotic coefficient and (D) activity coefficient of NaCl(aq) as a function of molality and temperature , as predicted from Clarke’s calculated virial matrix (surface plot) and reference data from the literature (Gibbard et al., 1974; Messikomer and Wood, 1975; Tanner and Lamb, 1978)

Reduced virial matrix and fit with experimental data

The problem with the virial matrix from Clarke and Glew (1985) is the high experimental cost associated with obtaining the virial coefficients, where measurements from more than five experimental techniques are required in order to calculate accurately all Virial coefficients up to order 6. In the present section, we will demonstrate that the use of a reduced virial matrix as shown in Equation (27) and Figure 3, which uses 20 virial coefficients instead of 32, results in an acceptable loss in accuracy for the estimation of osmotic and activity coefficients. As will be demonstrated in Section 4, the main interest of using the reduced virial matrix is the significant reduction in the experimental cost associated with the fitting procedure of the virial coefficients.

Figure 3

Graphical illustration of the fitting procedure for the reduced virial matrix

The reduced virial matrix is of the formwhere the coefficients are given by the Debye–Hückel theory as computed by Clarke and Glew (1980) and the coefficients are Pitzer coefficients fitted from experimental data of the excess heat capacity, enthalpy of mixing and osmotic coefficient at 25°C. As shown below, the reduced virial matrix allows a good estimation of activity coefficients, for three reference aqueous electrolytes: NaCl, KCl, and CaCl2.

The fitting datasets used for calculating the virial coefficients (shown in Table 1) were taken from Clarke and Glew (1985) for NaCl, from Parker (1965), Snipes et al. (1975) for KCl, and from Ananthaswamy and Atkinson (1985) for CaCl2.

Table 1

Coefficients of the reduced virial matrix fitted to available experimental data on osmotic coefficients, excess enthalpy and heat capacity from the literature

Solvent-specific	A_0	A_1	A_2	A_3	A_4	Ref.
	−116.8569	59.2284	0.772533	0.0113	1.3439 × 10−4	Clarke and Glew, 1980
	kg^12.mol^−12K	kg^12.mol^−12K	kg^12.mol^−12	kg^12.mol^−12K1	kg^12.mol^−12K2
Fit φ(25°C)	Q_0	B_0	C_0	D_0	E_0	Ref.
NaCl	−82.9	−22.51	−0.7836	0.09402	4.742e-3	Clarke and Glew, 1985
KCl	−73.02	−11.97	−0.9969	0.2107	−9.914e-3	Snipes et al., 1975
CaCl2	−497.4	−92.6	0.9724	−0.5565	0.05521	Ananthaswamy and Atkinson, 1985
	kg mol1K	kg mol1K	kg2mol2K	kg3mol3K	kg4mol4K
Fit ΔH^φ(25°C)	Q_1	B_1	C_1	D_1	E_1	Ref.
NaCl	36.7	81.84	−20.5	2.17	−0.1262	Clarke and Glew, 1985
KCl	80.02	61.8	−7.998	−0.1125	0.04836	Snipes et al., 1975
CaCl2	278.1	15.8	−6.714	−3.107	0.3062	Ananthaswamy and Atkinson, 1985
	kg mol1K	kg mol1K	kg2mol2K	kg3mol3K	kg4mol4K
Fit Cpφ(25°C)	Q_2	B_2	C_2	D_2	E_2	Ref.
NaCl	−0.4139	−1.723	0.4449	−0.05627	2.488e-3	Clarke and Glew, 1985
KCl	−2.028	−0.5265	−0.4449	0.2421	−0.03133	Tanner and Lamb, 1978
CaCl2	0.1441	−0.6466	−0.3304	0.07392	−3.602e-3	Ananthaswamy and Atkinson, 1985
	kg mol1	kg mol1	kg2mol2	kg3mol3	kg4mol4

The expressions at 25°C from (Equation 15), (Equation 16), (Equation 17) were used to obtain the values of the reduced virial coefficients shown on Table 1. The fits gave good results, with adjusted R-square >0.9996 for all electrolytes using MATLAB R2020b curve fitting tool. Once the virial coefficients were obtained, osmotic and activity coefficients for NaCl(aq), KCl(aq) and CaCl2(aq) could be calculated at different compositions and temperatures from Equation (15)-Equation (19), where the virial matrix was replaced with the reduced virial matrix introduced in Equation (27).

The obtained predicted values for the osmotic and activity coefficients for NaCl(aq), KCl(aq) and CaCl2(aq) are shown on Figures 4 and 5 (lines) together with the validation datasets (dots). The residuals (in %) between predicted and references values of the osmotic and activity coefficients for NaCl(aq), KCl(aq) and CaCl2(aq) are shown in C.

Figure 4

Reference literature values and calculated values for (A) the osmotic coefficient and (B) the activity coefficient of NaCl(aq) as a function of molality and temperature

Figure 5

Agreement between reference literature values and those calculated from the reduced virial matrix for (A) the osmotic coefficient and (B) the activity coefficient for KCl(aq), and for (C) the osmotic coefficient and (D) the activity coefficient for CaCl2(aq)

The validation dataset for NaCl(aq) was obtained by averaging reported values from Partanen and Partanen (2020), Gibbard et al. (1974) and Clarke and Glew (1985). Although the values reported by Partanen and Partanen (2020) were not in excellent agreement with the other reference values, a lower weight was given to this dataset.

It is apparent from Figure 4 that despite the reduction in the number of fitting coefficients, the agreement between predicted and reference values is very good on a range of molality spanning 0–5 mol/kg and a range of temperature spanning 0–60°C for NaCl(aq). The discrepancies between predicted values at different temperatures and reference literature values are given in C and show that the residuals are within [-0.43%, +0.49%] for the osmotic coefficient and [-0.92%, +1.05%] for the activity coefficient of NaCl(aq).

The validation set for KCl(aq) was taken from Snipes et al. (1975), Hostetler et al. (1967), Partanen (2016) and Hamer and Wu (1972), where, again, the dataset from Partanen et al. was given a lower weight. It is notable from Figure 5B that the calculated values are also in very good agreement with values reported in the literature, although, unfortunately, reference values for the activity coefficient of KCl(aq) at lower temperatures (0–25°C) could not be found. The discrepancies between predicted values at different temperatures and reference literature values are given in C and reveal that the residuals are within [-0.3%, +0.3%] for the osmotic coefficient and [-0.55%, +0.43%] for the activity coefficient of KCl(aq).

The validation set for CaCl2(aq) was taken from Ananthaswamy and Atkinson (1985). The agreement, although less satisfactory than for the previous 1:1 electrolytes, is still good: From C we can see that the residuals are within [−0.53%, +0.9%] for the osmotic coefficient and [−1%, +2%] for the activity coefficient. This is encouraging for the application of this method to aqueous organic 1:1 and 2:1 electrolytes such as the TEMPTMA - Methyl viologen system (Mourouga et al., 2022; Janoschka et al., 2016).

Measurement of activity coefficients

In the previous section, it was demonstrated that with the knowledge of the excess enthalpy , the excess heat capacity and the osmotic coefficient at a given temperature, the activity coefficient and osmotic coefficient could be calculated at any temperature with relatively good accuracy. In the present section, it will be demonstrated how differential scanning calorimetry (DSC) constitutes a rapid and efficient experimental method in the measurement of these parameters.

Differential scanning calorimetry

Differential scanning calorimetry is an experimental method where the temperature difference between an aluminum crucible containing the sample of interest (e.g., a 10 mg drop of electrolyte) and an empty, reference crucible is measured as a function of time. A linear temperature profile is imposed on the reference crucible (Figure 6A, red curve), so that when the sample holder reaches the freezing point, the exothermal freezing event results in a temperature increase in the sample holder (Figure 6A, blue curve). Similarly, on heating back to room temperature, the endothermal melting event results in a temperature decrease in the sample holder, which can be recorded and converted to units of mW/mg (Figure 6B).

Figure 6

Illustration of a NETZSCHE differential scanning calorimeter and (A) the temperature ramp imposed on the reference crucible and the temperature recorded on the sample holder and (B) the conversion to units of mW showing exothermal freezing and endothermal melting

It should be noted that the sample being unstirred during the experiments and the solvent being of high purity (mQ water), supercooling may occur during freezing events, so that the freezing temperature may be much lower than expected (pure water, for example, may freeze at −15°C). Supercooling is not observed on melting events, making them better suited for thermodynamic analyses.

On melting events, the location of the power peak is proportional to the melting temperature (in °C), the area under the peak to the enthalpy of melting (in J/g) and the difference between the baseline to the difference in sample heat capacity (in J/g/K) between solid and liquid state.

A method for measuring osmotic coefficients, that of freezing point depression (Lewis and Randall, 1923; Scatchard et al., 1932) allows to calculate the osmotic coefficient at the freezing point from the measurement of the freezing point of a mixture as a function of composition. Although DSC allows to measure the freezing point as well as the enthalpy and the heat capacity at the freezing point, it should be possible to derive , and from experiments, where is the freezing point of electrolyte solution. Following the reduced virial matrix method demonstrated in the previous sections, these measurements allow to calculate the activity and osmotic coefficients as a function of composition and temperature.

Illustration on a pure substance

Differential scanning calorimetry can be performed on a 10-mg drop of pure, mQ water to illustrate the experimental outputs of the method. Figure 7 shows screenshots of DSC results from experiments performed at LRCS, Amiens, France, where it should be noted that endothermal events are represented with a positive value on the axis. This is motivated by the supercooling observable on Figure 7A, where pure water freezes around −13°C, which makes freezing events unreliable, so that in the following graphs of Figure 7 only endothermal events are shown and the axis is shifted accordingly to manipulate positive quantities.

Figure 7

Graphs showing experimental DSC results on a 10 mg drop of mQ water. Endothermal events are represented with positive values

Figure 7B shows the measurement of the freezing point of pure water, which is defined as the intersect between the solid-state baseline (corresponding to the heat capacity of ice) and the tangent to the first inflection point of the peak (corresponding to the maximum variation of heat capacity during the phase transition). The intersect between these slopes defines the start of the phase transition, i.e., the melting point, where DSC exhibits a 0.67% accuracy. Figure 7C shows the measurement of the area under the power curve during melting, which should correspond to the enthalpy of melting for a pure substance. For pure water, the standard enthalpy of melting is reported as 333.25 J/g, so that DSC exhibits a 1.2% accuracy. Figure 7D shows the measurement of the difference between the baselines, which should correspond to the difference between the heat capacity of the solid and the liquid phase. For pure water, the heat capacity of ice is reported as 2.1 J/g/K and the heat capacity of water as 4.2 J/g/K, so that DSC exhibits a 9.5% accuracy. It should therefore be noted that DSC is not the most accurate experimental method, and that careful uncertainty propagation should be conducted over numerical calculations parametrized with DSC results.

Illustration on a mixture

Figure 8A shows the melting peaks of DSC experiments performed on aqueous sodium chloride solutions at different compositions. Compared to a pure substance, two melting peaks are visible: the left-most one corresponds to the eutectic point of NaCl(aq) (or, equivalently, to the melting point of NaCl(sat)). The second melting point corresponds to that of the ”free” water in the mixture and is therefore related to the melting point depression of NaCl(aq) as a function of composition.

Figure 8

Melting curves for NaCl(aq) obtained via DSC (A) at different compositions for a heating rate = 5 K/min and (B) at different heating rates for a composition of 2 mol/kg

Similarly to a pure substance, the area under the curve corresponds to the enthalpy of melting and the difference between the baselines to the difference in heat capacity between solid and liquid. Mixing between the saturated solution and the free water during the melting process allows to measure and from the DSC measurement, knowing the characteristics of the pure substances.

Limitations of the study

In the present study, it was demonstrated that activity and osmotic coefficients of aqueous mixtures of NaCl(aq), KCl(aq) and CaCl2(aq) could be derived with satisfying accuracy, provided literature data on their osmotic coefficient, enthalpy of mixing and excess heat capacity were available at a reference temperature. In order to apply these methods to novel electrolytes, an experimental method using differential scanning calorimetry (DSC) was demonstrated in principle. This method, however, may suffer from the relatively high uncertainties associated with DSC.

The main difficulty arises for the measurement of the melting point of aqueous mixtures, compared to that of pure substances: on the second peaks of Figure 8, the inflection point method cannot be used as the second melting peak does not exhibit an inflection point. Furthermore, Figure 8B shows the DSC curves at different heating rates. It is apparent that the heating rate has an impact on peak position and width, and therefore introduces uncertainties on the measurement of , and , due to thermal delays between the sample and the temperature sensor at the heating station.

In order to derive accurate values from the experiments, a deconvolution algorithm, allowing to reconstruct the theoretical curve at infinitely slow heating rate = 0 should be applied. This will be the focus of our next publication focused on presenting experimental deconvoluted differential scanning calorimetry results on reference electrolyte (NaCl, KCl, and CaCl2), as well as aqueous organic electrolytes for use in flow batteries (TEMPTMA chloride and methyl viologen dichloride).

Conclusion

In the present work, it was first shown that the field of flow battery modeling is increasingly questioning the validity of the dilute solution hypothesis, i.e., the fact that aqueous electrolytes behave ideally in the range of ionic strength suitable for flow battery operation. To correct this assumption, the most straightforward method is to use the concept of activity coefficient first introduced by Lewis (1907). However, very few authors provide estimated values of these activity coefficients for the aqueous electrolytes of interest in flow batteries, as a function of composition and, preferably, temperature.

We therefore demonstrated that the thermodynamic activity model developed by Pitzer (1973) and Clarke and Glew (1985) allows to calculate activity coefficients with excellent accuracy on a wide range of composition and temperature. The limitation, however, of the full Pitzer–Clarke method is the very high experimental cost associated with the calculation of all Pitzer coefficients.

A new method was therefore developed, namely the reduced virial matrix estimation, in order to reduce the experimental cost while keeping a satisfactory accuracy. The accuracy of the method was quantified with available reference literature data on activity and osmotic coefficients of common chloride electrolytes (namely NaCl(aq), KCl(aq) and CaCl2(aq)).

An experimental method using differential scanning calorimetry, based on the reduced virial matrix, was presented, and its challenges outlined. In order to apply the method to aqueous flow battery electrolytes, for which activity and osmotic coefficients have not yet been reported in the literature, a deconvolution algorithm will be developed in a subsequent publication. This combined theoretical and experimental method will allow for an accurate estimation of activity coefficients for parametrizing flow battery models, at a low experimental cost.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources, materials, data and code should be directed to and will be fulfilled by the lead contact, Gaël Mourouga (gaelmourouga@protonmail.com).

Materials availability

This study did not generate new unique reagents.

Data and code availability

The publicly available data analyzed in this paper has been compiled from the literature and converted to.xlsx format. Files are available on Mendeley data: https://data.mendeley.com/datasets/g4gmrwc9kg. The DOI is listed in the key resources table, and individual publications DOI are listed in the Reference section.

All original code has been deposited at Zenodo and is publicly available as of the date of publication. DOIs are listed in the key resources table. Executing the Matlab files, e.g. VirialMatrixNacl.m, requires some input data indicated in the first lines of the program. Note that these files should be in .txt format if the data was retrieved in .xslx format.

Both data and original code can be retrieved from the GitHub page https://github.com/Anon75014/Thermodynamics. The Matlab files should compile directly, provided the whole folder is downloaded at once.

Any additional information required to reanalyze the data reported in this paper is available from the lead contacton request.

Quantification and statistical analysis

Fitting of the reduced Virial matrix coefficients were conducted on MATLAB R2020b curve fitting tool (available on Zenodo) from selected experimental data from the literature (available on Mendeley Data). All fits yielded adjusted R-square values above 0.9996, up to 1 in most cases.

Gibbs free energy, chemical potential and activity coefficient

The Gibbs–Duhem equation states that at a specified temperature and pressure, the Gibbs free energy of an aqueous solution is proportional to the chemical potential of its individual components (Miller, 1963; Darken, 2002)where indexes + and - refer to cations and anions, respectively, and index w to water. The chemical potential of individual ions being a tricky quantity to define and manipulate (see Ebeling and Scherwinski (1983) and Malatesta (2010) for a lengthier discussion on the topic) we will conveniently introduce the chemical potential of soluteand the number of moles of soluteso thatwhich is Equation (12) in Section 2.3.

For a solute dissolved in water the Gibbs–Duhem equation reads

Kinetic theory of gases and ideal solutions

According to the kinetic theory of gases developped by Maxwell and Boltzmann in the 1870s, the chemical potential of non-interacting gas molecules (ideal gases) is proportional to the logarithm of their partial molar pressure where is the chemical potential of the gas in standard state, is the ideal gas constant, the absolute temperature and a reference partial pressure defined in the standard state. The same theory, applied to non-interacting ions in aqueous solutions (ideal aqueous solutions), would lead to a dependency of the chemical potential of solute in the logarithm of the molality (in moles of solute per kilogram of solvent)where represents a conventional reference molality of 1 mol/kg and the chemical potential for an ideal solution in standard state. The differential from of the Gibbs–Duhem equation Equation (A.5), at constant temperature and pressure, yieldsso that by differentiating Equation (A.7) and noting thatwhere is the molar mass of water (expressed in kg/mol), one can obtain the following expression for the chemical potential of water

The Gibbs free energy available in an ideal solution can therefore be calculated by inserting Equation (A.7) and Equation (A.10) in Equation (12) and subtracting the Gibbs free energy in standard state

Activity and chemical potential in non-ideal solutions

In the case of a non-ideal solution, the chemical potential is not directly proportional to the logarithm of the molality , but rather to the activity , as introduced by Lewis (1907)

The activity, as mentioned before, corrects the molality of solute by introducing an activity coefficient which accounts for the interactions between ions, as well as between the solute and the solvent. The activity is therefore defined from the molality aswhich leads to the following expression for the chemical potential of real solutions

By introducing the ideal chemical potential previously defined in Equation (A.7), it is possible to express the chemical potential of real solutions as a sum of an ideal term (corresponding to the hypothetical case of ideal solutions) and an excess term expressing the experimentally observed deviations in the behavior of real aqueous solutions

Similarly, the chemical potential of water can be expressed as the sum of an ideal term defined in Equation (A.10) and an excess term . The excess term, in turn, is proportional to the osmotic coefficient as

Gibbs free energy of real solutions

Defining the Gibbs free energy of the solution as the sum of an ideal component , corresponding to the potential energy which can be extracted from a hypothetical ideal solution, and an excess component , corresponding to the deviations observed from the hypothetical ideal case, we therefore obtainwhereand where the partial molal Gibbs free energy is simply defined by normalizing the excess term with the number of moles of solutewhich yields Equation (14) introduced in Section 2.3

Fundamentally, the excess partial molal Gibbs free energy captures all deviations from the hypothetical case of ideal solutions (as predicted by the kinetic theory of gases) observed in colligative properties of real aqueous solutions.

MatrixForm of the virial and Taylor expansions

The full matrix form of (Equation 21), (Equation 22), (Equation 23), (Equation 24) can be explicited as:

Residuals between literature references and predicted values using the reduced virial matrix

Figure 9: Residuals between reference literature values (Gibbard et al., 1974; Clarke and Glew, 1985; Partanen and Partanen, 2020) and those calculated from the reduced virial matrix for (a) the osmotic coefficient (ϕ) and (b) the activity coefficient (γ) of NaCl(aq) as a function of molality and temperature.

Figure 10: Residuals between reference literature values literature (Hostetler et al., 1967; Hamer and Wu, 1972; Snipes et al., 1975; Platford, 2002; Amado G and Blanco, 2004; Partanen, 2016) and those calculated from the reduced virial matrix for (a) the osmotic coefficient (ϕ) and (b) the activity coefficient (γ) for KCl(aq).

Figure 11: Residuals between reference literature values (Ananthaswamy and Atkinson, 1985) and those calculated from the reduced virial matrix for (a) the osmotic coefficient (ϕ) and (b) the activity coefficient (γ) for CaCl2(aq).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Software and algorithms
MATLAB code	Database: Zenodo	https://doi.org/10.5281/zenodo.6926800
Deposited data
Experimental data compiled from lit.	Database: Mendeley Data	https://doi.org/10.17632/g4gmrwc9kg.1