Thermal conductivity of aqueous solutions of reline, ethaline, and glyceline deep eutectic solvents; a molecular dynamics simulation study

Abstract

Accurate knowledge and control of thermal conductivities is central for the efficient design of heat storage and transfer devices working with deep eutectic solvents (DESs). The addition of water is a straightforward and cost-efficient way of tuning many properties of DESs. In this work, the thermal conductivities of aqueous solutions of reline, ethaline, and glyceline are reported for the first time. The non-equilibrium molecular dynamics Müller-Plathe (MP) method was used, along with the well-established GAFF and SPC/E force fields for DESs and water, respectively. We show that thermal conductivities of neat DESs are in excellent agreement with available experimental data. The addition of 25 wt% water results in nearly 2 times higher thermal conductivities in all DESs. A further increase in the fraction of water to 75 wt%, causes an increase in the thermal conductivities of DESs ca. 3 times. This behaviour is mainly due to the change in the microscopic structure of the DESs (i.e. hydrogen bonding) upon the addition of water. Our simulations reveal that thermal conductivities of aqueous DESs do not significantly depend on temperature. We also show that thermal conductivities strongly depend on system-size. System-sizes bigger larger than ca. 5 nm should be used.

GRAPHICAL ABSTRACT

Keywords:

• Thermal conductivity
• deep eutectic solvents
• molecular dynamics simulation
• aqueous solutions
• temperature effects

1. Introduction

Sustainable, efficient, and cost-effective heat transfer fluids are essential for a wide variety of industrial and domestic applications such as solar power generation [1, 2], heat storage in power plants [3], high-temperature processing of plastics [4], refrigeration systems [5], and cryogenic gas processing [6, 7]. The conventional organic solvents used in these processes, e.g. benzene, chloroform and toluene, often exhibit non-desirable properties such as high-toxicity and/or high-volatility [8, 9]. To this purpose, the demand for non-hazardous solvents with tailored properties is great. Deep Eutectic Solvents have emerged as a new generation of environmentally-friendly solvents, which can be tuned to exhibit enhanced physical properties [10–12].

DESs are mixtures of a Hydrogen Bond Acceptor (HBA) (e.g. choline chloride, tetrabutylammonium bromide) and a Hydrogen Bond Donor (HBD) (e.g. urea, malonic acid, glycerol) in specific ratios. DESs are liquids having melting points considerably lower than the ones of the building blocks used for their synthesis [13–16]. Properties such us the low volatility, low vapour pressure, large liquid range, large heat capacity, high thermal stability, and non-flammability, make DESs very promising working fluids in high-temperature applications [13, 17–19]. Depending on the nature of the building blocks used, DESs can be biodegradable, bio-compatible, and non-toxic [20–22]. For this reason, DESs can be also used in low-temperature refrigeration systems [23], and in the pharmaceutical industry [24]. The relatively low-cost of the starting materials and the easy synthetic pathway may enable large-scale industrial production of DESs [21].

The accurate prediction of the thermal conductivities of DESs is crucial for the safe and efficient design of heat transfer processes and devices working with DESs (e.g. chemical reactors, heat exchangers, thermal storage units) [25–27]. Thermal conductivity is a measure for the amount of heat that is transferred through a material due to a temperature gradient [25, 26]. Despite being a key property, only a very limited number of studies in literature have focussed on the thermal conductivity of DESs. Recently, the thermal conductivities of choline chloride (ChCl)/urea [26], ChCl/glycerol [28, 29], and ChCl/ethylene glycol [30] have been measured experimentally at different temperatures. It was shown that ChCl-based DESs have relatively low thermal conductivities that do not strongly depend on temperature. The low thermal conductivities, combined with the high viscosities, present a significant advantage for many practical applications of DESs, for example, as a thermal barrier (i.e. insulating materials) in devices that fail due to overheating [31]. A low thermal conductivity is also considered to be a limitation in the efficiency of heat transfer devices such as heat sinks in electronics [32]. To this purpose, tailor-designed DESs having thermal conductivities that can satisfy certain applications is pivotal. A relatively easy and cost-effective way to tune the thermodynamic and transport properties of DESs is by adding controlled amounts of water [33–35]. As shown by Shah and Mjalli [34], Celebi et al. [35], Fetisov et al. [36], and Kumari et al. [37], the addition of water significantly changes the molecular structure of DESs (i.e. hydrogen bond network). As a result, a number of properties of DESs such as the densities [6, 33, 34, 38–41], viscosities [6, 33, 34, 38–40, 42–44], diffusivities [6, 34, 35, 44], ionic conductivities [34, 35, 40], speed of sound [34, 42], and solubilities [16, 45, 46] are vastly affected, depending on the water content in the mixture.

To the best of our knowledge, no experimental or computational work in literature has focussed on the thermal conductivities of aqueous DESs solutions. To this end, we present new data for the thermal conductivities of aqueous reline, ethaline, and glyceline solutions computed using non-equilibrium molecular dynamics (NEMD) for a wide range of mass fractions of water. The temperature-dependence of thermal conductivities of neat and aqueous glyceline solutions for temperatures ranging from 283.15 to 363.15 K is also shown. A detailed discussion of the effects of the system size, and the length of the swap integral in the MP method on the computation of thermal conductivities is also provided. The remainder of the paper is organised as follows: In Section 2, the description of the used force fields, simulation details, and the MP method are presented. Our main findings are presented in Section 3. Finally, conclusions are discussed in Section 4.

2. Model and method

2.1. Force fields

The generalised amber force field (GAFF) was used to model reline, ethaline, and glyceline [47]. The partial atomic charges were derived via the Restrained Electrostatic Potential (RESP) method [48–50], using the Hartree-Fock HF/6.31G* level of theory [49]. The partial charges of ChCl in reline were uniformly scaled by a factor of 0.8 [49, 51]. The partial charges of ChCl in ethaline and glyceline were uniformly scaled by a factor of 0.9 [51]. As discussed in detail in the studies of Perkins et al. [49, 51] and Liu et al. [18], charge scaling in ILs and DESs is essential because electrostatic interactions are usually overpredicted. Recently, Blazquez et al. [52] showed that the scaling of charges of even simple electrolytes such as NaCl can lead to an improved prediction of the salting out effect of methane in water. In the past few years, various studies have shown that simulations using the GAFF force field combined with scaled charges yield relatively accurate thermodynamic and transport properties of neat reline, ethaline, and glyceline [6, 16, 34, 35, 49, 51, 53].

Water was modelled using the rigid three-site SPC/E force field [54]. SPC/E has been shown to be accurate in predicting the transport properties of pure water [55, 56], and mixtures of water with gases [57] for a wide range of temperatures. The SPC/E force field has been also widely used to model aqueous NaCl solutions [58–63], and aqueous solutions of ionic liquids (ILs) and DESs [6, 34, 64]. In the case of DESs, the combination of GAFF and SPC/E force fields have been shown to yield good agreement with experiments for several thermodynamic and transport properties of aqueous reline, ethaline, and glyceline solutions such as viscosities, diffusivities, ionic conductivities, and thermal expansivities [6, 34, 35]. For all simulations of argon (i.e. for the system size dependence study of thermal conductivity) the parameters of Ghorbanian et al. [65] are used. All force field parameters and charges used in our simulations are listed in Tables S1–S13 of the Supplementary Material.

To further validate the performance of the chosen force fields, we performed equilibrium molecular dynamics (EMD) simulations to compute the densities of aqueous DES solutions. These values were compared to the available experimental data [33, 38, 39]. As can be seen in Figure S4 of the Supplementary Material, the computed densities are in good agreement with the experimental values (deviations are at most 1.5%). It is important to note that the scope of this work is not to exhaustively evaluate the performance of force field combinations, but to investigate the underlying molecular phenomena affecting the thermal conductivity of aqueous DES solutions, and provide new reliable data.

2.2. Molecular dynamics simulation details

MD simulations were carried out in a cubic simulation box, with periodic boundary conditions applied in all directions. Initial configurations of aqueous reline, ethaline, and glyceline solutions were generated using the Packmol software [66]. For each DES, eutectic compositions with a 1:2 molar ratios were used, in which 1 refers to the HBA, and 2 refers to the HBD [67]. To study the effect of water content on the thermal conductivity of the aqueous DES solutions, the mass fraction of water (${\omega }_{\mathrm{w}}$) was varied according to: (1) ${\omega }_{\mathrm{w}}=\frac{N_{\mathrm{w}}{M}_{\mathrm{w}}}{N_{\mathrm{ChCl}}{M}_{\mathrm{ChCl}}+N_{\mathrm{HBD}}{M}_{\mathrm{HBD}}+N_{\mathrm{w}}{M}_{\mathrm{w}}}$(1) where $M_i$ and $N_i$ are the molecular weight and the number of molecules of species i (i.e. ChCl, HBD, and water), respectively. In this study, ChCl is considered to be a single component. In all simulations, the system sizes correspond to simulation boxes equal to ca. 5.5 nm. This was determined by systematically increasing the simulation box lengths until no significant variations in the computed thermal conductivities of DESs are present. A detailed analysis of the system size dependence of thermal conductivities is provided in Section 3.2. A list of all simulated systems can be found in Table S14 of the Supplementary Material.

All MD simulations were performed using the open-source Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS), version released on August, 2018 [68]. For integrating Newton's equations of motion, the Verlet algorithm was used with a time step of 1 fs. Long-range electrostatic interactions between charged species were handled using the particle-particle, particle-mesh (pppm) solver with a root-mean precision of ${10}^{-5}$ [69]. The cut-off distance was set to 12 Å, for both Lennard-Jones (LJ) and the real-space part of electrostatic interactions. The LJ parameters for the interaction between dissimilar atoms were obtained using the Lorentz–Berthelot combining rules [70]. For water molecules, bond lengths and bond-bending angles were kept fixed using the SHAKE algorithm [71].

The simulation scheme used in all systems is described in the following steps:

1. Initial configurations were energy minimised for 10,000 steps using the conjugate gradient method.

2. EMD simulations of 1 ns in the isothermal-isobaric (NPT) ensemble were performed for the equilibration of the systems. The temperature and pressure were kept constant using the Nosé-Hoover thermostat and barostat with coupling constants of 0.1 ps and 1 ps, respectively.

3. NPT runs of 4 ns are performed for computing the average volumes and densities. The densities of all simulated aqueous DES solutions are listed in Table S15 of the Supplementary Material along with the respective experimental values.

4. Starting from the average densities obtained in the previous step, each system was further equilibrated for 1 ns in the canonical (NVT) ensemble. The temperature was kept constant using the Nosé-Hoover thermostat with a coupling constant of 0.1 ps.

5. NEMD production runs of 50 ns in the NVT ensemble were performed for computing the thermal conductivities using the MP method [72]. The MP method is discussed in detail in the following section. The standard deviation of the computed thermal conductivities are obtained from 6 independent simulations, each one starting from a different initial configuration.

2.2.1. Müller-Plathe (MP) method for computing thermal conductivities

Thermal conductivities can be computed by NEMD and EMD simulations. In NEMD, the response of the system to an externally induced heat flux, or an imposed temperature gradient, yields the thermal conductivity [72–77]. In EMD, two equivalent methods can be used, i.e. the Green-Kubo (GK) [73, 78–81] or Einsten relations [31, 82–84]. It is important to note that both EMD methods require the heat flux vector J for computing thermal conductivities [73]. As shown in the recent studies by Surbyls et al. [85] and Boone et al. [86], the computation of J in molecular simulation packages such as LAMMPS can be erroneous for molecular systems having angles, torsions, dihedrals, and improper dihedrals, producing unphysical values. In fact, our preliminary thermal conductivity computations for reline using the Einstein method (OCTP plugin in LAMMPS [84]) yielded values which were almost an order of magnitude higher than the respective experiments. Therefore, in this study, the MP method was used. The MP method [72] is based on NEMD simulations and has been extensively used for computing thermal conductivities of fluid models, e.g. Lennard-Jones (LJ) particles [72, 87], as well of a wide variety of molecular fluids and liquid mixtures, e.g. water, hydrocarbons and ionic liquids (ILs) [18, 67, 88, 89]. A brief description of the MP method for computing thermal conductivity is as follows: A heat flux is induced to the system by exchanging kinetic energies between the hottest and coldest particles in different regions of the simulation box. Energies are exchanged with a swap integral W (in units of time). As we discuss in detail in the Section 3.1, the choice of W governs the magnitude of the heat flux (and of the resulting temperature gradient). The time-averaged heat flux $⟨J_z(t)⟩$ is described by: (2) $⟨J_z(t)⟩=\frac{\sum _{\mathrm{transfers}}\frac{m}{2}(v_{\mathrm{hot}}^2-v_{\mathrm{cold}}^2)}{2t{L}_{\mathrm{x}}{L}_{\mathrm{y}}}$(2) where the sum is taken over all energy transfers during the simulation time t. $v_{\mathrm{cold}}$ and $v_{\mathrm{hot}}$ denote the velocities of the cold and hot atoms. The mass of the exchanged atoms are represented by m. The product $L_x{L}_y$ is the cross-sectional area in which heat transport occurs and z is the direction of the temperature gradient. The factor 2 in the denominator is due to the periodicity, indicating that energy can be transferred in two directions. The velocity exchange increases the local temperature of the hot slab while decreasing the local temperature of the cold slab. Since the total energy is conserved, this velocity exchange eventually creates a temperature gradient $\mathrm{d}T/\mathrm{d}z$ in the system. The temperature gradient is computed by linear regression of local temperatures in the flow direction for each half of the simulation box. The thermal conductivity (λ) is computed by the ratio of the imposed heat flux to the resulting temperature: (3) $\lambda =-\frac{⟨J_z(t)⟩}{⟨\mathrm{d}T/\mathrm{d}z⟩}$(3) where the brackets $⟨\cdots ⟩$ denote a time average.

3. Results and discussion

In this section, our findings for various aspects that can affect the thermal conductivity of aqueous DESs solutions are presented. We discuss the effect of the swap integral in the MP method, system size effects, the effect of water content, and the temperature dependence of thermal conductivity.

3.1. The effect of the swap integral in the MP method

In the MP method, the generated heat flux (and thus the resulting temperature gradient) strongly depends on the choice of the swap interval W [72]. To investigate the effect of W on the computation of thermal conductivities of aqueous DESs solutions, we performed simulations for different values of W . In Figure 1(a), the resulting temperature profiles as a function of distance in the z-direction (i.e. the direction of the temperature gradient) are shown for W = 20, 40, 80, 160 and 320 fs for neat reline. The temperature gradient decreases as W increases. For the lowest value of W tested (i.e. 320 fs) a temperature difference of 10 K between the hot and cold regions is obtained. In contrast, for low values of W (i.e. $W=25\mathrm{fs}$), the temperature difference between the hot and cold regions equals ca. 100 K. For the other values of W , the temperature difference between hot and cold regions lie between these two values. As shown in the studies by Müller-Plathe [72] and Sirk et al. [88], sufficiently high values of W are essential to generate large enough heat fluxes (and the resulting temperature gradients). The value of W should be low enough to ensure that linear response theory still holds (i.e. linear temperature profiles are obtained). Consequently, a proper selection of W is important to minimise the error in the computed thermal conductivities due to the linear regression of the temperature profiles. In Figure 1(b), we show that the computed thermal conductivities for different values of W are practically identical within the error bars. No significant differences in the computation time of the thermal conductivities were observed by using the different values of W . The smallest error (i.e. 2%) in the computed thermal conductivity is observed for $W=80\mathrm{fs}$. For this reason, $W=80\mathrm{fs}$ is used in all the remainder simulations in this work.

Figure 1. The effect of swap integral on the (a) resulting temperature gradient, and (b) computed thermal conductivities of neat reline solutions at 303.15 K and 1 atm.

The presence of water in the DES solution also effects the magnitude of the computed temperature gradient. In Figure 2, the temperature gradient of aqueous reline solutions at 303.15 K and 1 atm are shown as a function of distance in z-direction for different water contents using $W=80\mathrm{fs}$. It is clear that the temperature gradient decreases the water content in the solution increases. The temperature difference between the hot and the cold regions for the neat reline and neat water are 45 K and 15 K, respectively. For the rest of the mass fractions of water, the temperature difference between the hot and cold regions are between these two values. These profiles are adequate for linear regression, and thus, a precise computation of thermal conductivity. This further justifies the choice of $W=80\mathrm{fs}$.

Figure 2. Temperature gradients of aqueous reline solutions at 303.15 K and 1 atm as a function of distance in z-direction for mass fractions of water ranging from 0 to 100%. $W=80\mathrm{fs}$ is used in all simulations.

3.2. The effect of the system size on thermal conductivities of aqueous DESs

In molecular simulations, the computation of various transport and thermodynamic properties such as the self- and collective diffusivities [83, 90–93], Kirkwood-Buff integrals [94–97], activity coefficients [98], and ionic conductivities [99] strongly depends on the system size. This is mainly due to the use of periodic boundary conditions, and the long-range nature of hydrodynamic and electrostatic interactions [100, 101]. To mitigate these finite-size effects, the use of relatively large system sizes or the use of analytic corrections is crucial. Thermal conductivities computed with the MP method also exhibit finite-size effects [88, 102]. Sellan et al. [102] showed that these effects become significant when the size of the simulation box is smaller than the mean-free path of the phonons.

In Figure 3, the effect of system size on the computed thermal conductivities of liquid argon at 86 K and 1 atm (which corresponds to an average density of $1403{\mathrm{kg}/\mathrm{m}}^3$), and neat reline at 303.15 K and 1 atm (which corresponds to an average density of $1212{\mathrm{kg}/\mathrm{m}}^3$) are shown. Argon was used as a test-case because it is very computationally efficient (i.e. modelled as a single LJ interaction site) allowing for the investigation of big system sizes. The largest system size of argon in this work is 571,787 molecules, which corresponds to a simulation box length ($L_c$) of ca. 30 nm. Performing simulations with similar box lengths for reline is very computationally intensive. The largest tested system size for reline is 1575 molecules, which corresponds to a simulation box of ca. 8 nm. For both argon and reline, the computed thermal conductivities initially increase with the system size until converging to a plateau at a system size of ca. 5 nm, which corresponds to systems containing more than 2000 and 400 argon and reline molecules, respectively. As shown in Figure 3, the computed thermal conductivities (for system sizes at which the plateau is reached) of both argon and reline are in good agreement with experimental results, showing a deviation of ca. 5%. Based on this investigation, all simulations for the computation of thermal conductivities are performed using system sizes corresponding to simulation box lengths of 5.5 nm. The number of molecules used in each simulation is listed in Table S16 in the Supplementary Material.

Figure 3. System size dependence of computed thermal conductivities of (a) liquid argon at 86 K and 1 atm ($\rho =1403{\mathrm{kg}/\mathrm{m}}^3$), and (b) neat reline at 303.15 K and 1 atm ($\rho =1212{\mathrm{kg}/\mathrm{m}}^3$). $N_{\mathrm{Ar}}$ and $N_{\mathrm{R}}$ are the number of argon and reline molecules, respectively. The red dotted lines refer to experimental data for argon [119] and reline [26]. The dashed lines connecting the symbols are to guide the eye.

3.3. The effect of the water content on thermal conductivities of aqueous DESs

In Figure 4, the computed thermal conductivities of aqueous reline, ethaline, and glyceline solutions at 303.15 K are shown as a function of the mass fraction of water. To the best of our knowledge, no experimental measurements of the thermal conductivity of aqueous DES solutions are available, while a few are reported for the neat DESs [26, 28–30]. As shown in Figure 4, computed thermal conductivities of neat reline, ethaline, and glyceline, are equal to $0.244\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$, $0.228\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$, and $0.237\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$, respectively. These results are in good agreement with the available experiments, i.e. the thermal conductivities of neat reline at 298 K and neat glyceline at 305 K were reported to be $0.245\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$ and $0.232\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$, respectively [26, 29]. Interestingly, the thermal conductivities of many common ILs, such as [EMIM][OAc], [BBIM][NTf${}_2$] [C${}_2$mim][EtSO${}_4$], and [C${}_4$mim][Tf${}_2$N] are in the range of 0.11 and $0.23\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$ [18, 89, 103, 104]. Nevertheless, the thermal conductivities of DESs and ILs are much lower compared to the ones of aqueous solutions of simpler electrolytes such as NaCl, LiBr, and CoCl${}_2$ which are in the range $0.5--0.7\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$ [105–107].

Figure 4. The thermal conductivities of aqueous reline, ethaline, and glyceline solutions at 303.15 K and 1 atm as a function of the mass fraction of water. The symbol × represents the computed thermal conductivity of the quaternary mixture of reline+etaline+glyceline+water in which each component has a mass fraction of 25 wt%. BlUe, green and red lines (dashed and solid) represent the fit using the Jamieson correlation (Equation (4(4) ${\lambda }_m={\omega }_1{\lambda }_1+{\omega }_2{\lambda }_2-\alpha ({\lambda }_2-{\lambda }_1)[1-(\sqrt{{\omega }_2})]{\omega }_2$(4) )) for reline, glyceline and ethaline, respectively [25]. The uncertainties in the computed thermal conductivities are in the range of 1–3%.

As shown in Figure 4, the thermal conductivities of the aqueous DES solutions increase with the increase in the water content. This can be explained by the changes in the micro-structure of the DESs in the presence of water. As shown in previous studies [6, 34–37], with the addition of water, hydrogen bonds (HB) between HBDs and anions, and between anions are significantly reduced while new HBs between the water molecules and the ions are formed. Moreover, water reduces the viscosity of the mixture, and thus, the self-diffusivities of all species increase along with the thermal conductivities [108, 109]. More details on the variation of the HBs and self-diffusivities as a function of the water content can be found in the studies by Baz et al. [6] and Celebi et al. [35].

In this study, our computed thermal conductivity of the SPC/E water model is equal to $0.853\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$ (see Figure 4). This value is in agreement with earlier results for the same force field computed either using the MP [67, 88, 110, 111] or the Green-Kubo method [88, 112, 113]. The thermal conductivity of the rigid SPC/E water model at 303.15 K deviates from the experimentally measured one (i.e. $0.615\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$) by ca. 35%. Therefore, given the good performance of the GAFF force field in predicting the thermal conductivity of neat DESs, and the lower accuracy of the SPC/E force field for water, it is reasonable to expect that the computed thermal conductivities of the aqueous DES solutions will deviate from the actual thermal conductivities by not more than 35%, depending on the water content in the mixture. It is important to note that such overestimations of the thermal conductivity of water were also reported by many previous MD studies for various water force fields (i.e. 3-site, 4-site, flexible, and rigid) [67, 88, 110–113]. In these studies it is shown that none of the considered water force fields can accurately reproduce the experimentally measured thermal conductivity. Despite this deviation, it is important to note that our findings for the dependence of thermal conductivity of the aqueous DES solutions on the water content are in qualitative agreement with prior experimental and simulation studies on aqueous ILs solutions. Ge et al. [104] performed experiments to determine the effect of water content on the thermal conductivities of aqeuous solutions of [C${}_2$mim][EtSO${}_4$] and [C${}_2$mim][OTf] at 293 K. Kelkar et al. [114] performed MD simulations of aqueous [C${}_2$mim][EtSO${}_4$] solutions. In both studies, an increase in the thermal conductivities was obtained with increasing water content.

As can be seen in Figure 4, the computed thermal conductivities of the aqueous reline solutions are slightly, but consistently, higher than the respective ones of the aqueous glyceline and ethaline solutions for the whole range of water mass fractions [35]. These small differences can be mainly explained by the different H-bonding behaviour of the three DESs due to the variation in the HBDs, i.e. urea in reline, ethylene glycol in ethaline and glycerol in glyceline [6, 35]. Urea is a strong HBD, which has been experimentally shown to be associated with high thermal conductivities. Gautham and Seth [26] measured the thermal conductivities of ChCl/urea and ChCl/thio-urea and showed that ChCl/urea has higher thermal conductivity due to the stronger H-bonding nature of urea (i.e. the electronegativity of oxygen in urea is higher than sulfur in thio-urea). Our results also show that the thermal conductivities of aqueous glyceline solutions are slightly higher than those of ethaline. This is consistent with the experimentally measured thermal conductivities of the neat glycerol and etyhlene glycol. The measured thermal conductivity of glycerol ($0.292\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$) is slightly higher than the one of ethylene glycol ($0.256\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$) at 298.15 K [115] In this study, we also simulated a quaternary mixture of reline-ethaline-glyceline-water with equal mass fractions. We aimed to diversify HB networks using three different DESs in an aqueous solution, and investigate how this affects the thermal condutivities. The computed thermal conductivitiy of the quaternary mixture is equal to $0.413\mathrm{W}\cdot {\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$. As can be clearly seen in Figure 4, this value is nearly identical to the thermal conductivity of glyceline solution having mass fraction of water of 25 wt%, which lies between the thermal conductivities of reline and ethaline solutions for the same water content. This finding clearly shows that the thermal conductivities of DESs can be accurately tuned by adjusting the compositions in their mixtures. All raw data and the corresponding uncertainties are listed in Table S17 of the Supplementary Material.

Several empirical models such as the Flippov equation [25], Jamieson correlation [116], Baroncini correlation [25, 116], and Rowley method [117] have been proposed for predicting thermal conductivities of liquid mixtures from pure component data. Jamieson correlation [116] was validated using the thermal conductivities of a wide range of aqueous and non-aqueous binary mixtures for various compositions. These mixtures involved ethers, ketones, alcohols, aromatics, weak acids, water, saccharides, milk, butter, whey, and fruit juices. Jamieson correlation is as follows: (4) ${\lambda }_m={\omega }_1{\lambda }_1+{\omega }_2{\lambda }_2-\alpha ({\lambda }_2-{\lambda }_1)[1-(\sqrt{{\omega }_2})]{\omega }_2$(4) where ${\lambda }_m$ is the thermal conductivity of the binary mixture of components 1 and 2, having thermal conductivities ${\lambda }_1$ and ${\lambda }_2$, respectively. The two components are sorted so that ${\lambda }_2>{\lambda }_1$. $w_1$ and $w_2$ are the mass fractions of components 1 and 2, respectively. α is an empirical parameter that is temperature-independent and can be obtained by fitting to the available experimental data. According to Jamieson et al. [116], the thermal conductivity of any binary mixture (aqueous or non-aqueous) can be predicted with a minimum accuracy of 7% by Equation (4(4) ${\lambda }_m={\omega }_1{\lambda }_1+{\omega }_2{\lambda }_2-\alpha ({\lambda }_2-{\lambda }_1)[1-(\sqrt{{\omega }_2})]{\omega }_2$(4) ) with $\alpha =1$ (if no experimental data are available to fit).

Although Jamieson correlation was not developed considering aqueous electrolyte solutions, it has been used earlier to predict the thermal conductivities of such systems. Ge et al. [104], validated Equation (4(4) ${\lambda }_m={\omega }_1{\lambda }_1+{\omega }_2{\lambda }_2-\alpha ({\lambda }_2-{\lambda }_1)[1-(\sqrt{{\omega }_2})]{\omega }_2$(4) ) using the thermal conductivities of aqueous solutions of [C${}_2$mim][EtSO${}_4$] and [C${}_2$mim][OTf] ILs. These authors showed that by using fitted α parameters equal to 0.4211 and 0.7043, respectively, Equation (4(4) ${\lambda }_m={\omega }_1{\lambda }_1+{\omega }_2{\lambda }_2-\alpha ({\lambda }_2-{\lambda }_1)[1-(\sqrt{{\omega }_2})]{\omega }_2$(4) ) is able to predict the thermal conductivities very accurately. Fitting the data by Ge et al. [104] using the general rule proposed by Jamieson et al. [116], i.e. $\alpha =1$, yields thermal conductivities with accuracy in the range of 1.7–14.8%. Excluding the mass fractions of water equal to 20 wt% and 50 wt%, the predictions of thermal conductivities show an accuracy in the range of 1.7–6.6%. This clearly shows that Jamieson correlation can be safely used to predict the thermal conductivities of aqueous ILs solutions, even by using a fixed value of $\alpha =1$.

As shown in Figure 4, the computed data for the aqueous DESs solutions are fitted using the Jamieson correlation (Equation (4(4) ${\lambda }_m={\omega }_1{\lambda }_1+{\omega }_2{\lambda }_2-\alpha ({\lambda }_2-{\lambda }_1)[1-(\sqrt{{\omega }_2})]{\omega }_2$(4) )). Values of fitted α parameters equal to −0.61, −0.22 and 0, for reline, glyceline and ethaline solutions, respectively, yield a very good fit of the Jamieson correlation to the computed thermal conductivities. The thermal conductivity of aqueous ethaline is an almost linear function of the mass fraction of water, and thus, $\alpha =0$ in this case. In Figure 4, we also fit the computed data using $\alpha =1$. The use of $\alpha =1$ does not yield accurate predictions for the thermal conductivities of the aqueous DESs solutions, showing deviations from the computed data in the range of 8.8–30.4%. Interestingly, if the fitted α parameters are used, the Jamieson correlation is a concave function, while the use of $\alpha =1$ produces a convex function of the thermal conductivities with the mass fraction of water. These findings clearly show that the general rule of thumb proposed by Jamieson et al. (i.e. $\alpha =1$) does not hold for the prediction of thermal conductivities of aqueous DESs solutions, and hence new models may be needed.

3.4. The effect of the temperature on thermal conductivities of aqueous DESs

Figure 5 shows the computed thermal conductivities of aqueous glyceline solutions for temperatures in the range of 283.15–363.15 K, and mass fractions of water from 0 (i.e. neat glyceline) to 100 wt% (i.e. neat water). For the whole range of mass fractions of water, it is evident that the thermal conductivities do not depend on the temperature. This can be mainly explained by the negligible change in the microstructure of aqueous DESs (and neat water) with the increase in temperature. As recently shown by Celebi et al. [35], the radial distribution functions and the HB networks between the components of aqueous reline and ethaline solutions do not significantly change as a function of temperature. These findings are also in line with several earlier experimental measurements for different DESs [26, 29, 30], and MD studies and experiments for many common ILs [18, 89, 103, 104, 118]. The thermal conductivities of ILs showed a small decrease with the increase in temperature [18, 89, 103, 104, 118]. In contrast, the experimental study by Singh et al. [28] revealed that a temperature increase of ${80}^{\circ }\mathrm{C}$ results in a more than 10% increase in the thermal conductivity of neat glyceline. As shown in Figure 5, for the same system a temperature increase from 283.15 to 363.15 K results in 4% increase in the thermal conductivity. Note that this difference is within the uncertainty range (see Table S18 of the Supplementary Material).

Figure 5. The effect of temperature on the computed thermal conductivities of aqueous glyceline solutions with water mass fractions ranging from 0 to 100%. The uncertainties in the computed thermal conductivities are in the range of 1–3%. The dotted lines connecting the symbols are to guide the eye.

4. Conclusions

We performed molecular dynamics simulations using the Müller-Plathe method [72] to compute thermal conductivities of aqueous solutions of reline, ethaline, and glyceline DESs at 303.15 K, and mass fractions of water ranging from 0 to 100 wt%. To the best of our knoweledge, no available experimental or simulation data are available for the aqueous DESs mixtures. The GAFF and SPC/E force fields were used for modelling DESs and water, respectively. Our findings for neat DESs show a good agreement with available experimental data, deviating by a maximum of 3%. With the addition of 25 wt% water, thermal conductivities of all DESs increase by almost a factor of 2 (i.e. 1.9 for reline, 1.8 for glyceline, and 1.7 for ethaline). For 75 wt% water, the thermal conductivities increase ca. 3 times. This behaviour is mainly related to the depletion of hydrogen bonds between the components of DESs, and the formation of new ones, between the water molecules and DESs. The differences between the computed thermal conductivities of the three DESs are mainly due to the hydrogen bonding ability of urea, ethylene glycol, and glycerol. Our simulations revealed that the influence of temperature on the thermal conductivities of the aqueous solutions of DESs is negligible. For the aqueous glyceline solutions, the computed thermal conductivities at 283.15 K are ca. 4% lower than the ones at 363.15 K for the whole range of mass fractions of water. We also investigated the dependence of the computed thermal conductivities on the simulation system size using argon and reline as test cases. The thermal conductivities of both systems showed an initial increase with the system size until reaching a constant value for system sizes corresponding to simulations box lengths of ca. 5 nm. The effect of the swap integral in the Müller-Plathe method was also studied. Although the swap integral governs the magnitude of the heat flux, we show that all values in the range of 20–320 fs yield practically identical thermal conductivities. Nevertheless, the small and large values cause relatively large errors in the computed thermal conductivities. Therefore, the use of an intermediate swap integral, e.g. 80 fs, is advised for molecular simulations of aqueous DESs solutions.

Acknowledgements

This study was funded by Nederlandse Organisatie voor Wetenschappelijk Onderzoek within the framework of NWO Exacte Wetenschappen (Physical Sciences) for the use of high-performance computing facilities. TJHV acknowledges NWO-CW (Chemical Sciences) for a VICI grant.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This study was funded by Nederlandse Organisatie voor Wetenschappelijk Onderzoek within the framework of NWO Exacte Wetenschappen (Physical Sciences) for the use of high-performance computing facilities. TJHV acknowledges NWO-CW (Chemical Sciences) for a VICI grant [NWO-CW,NWO-EW].