Thermal Transport Anomalies of Electrolyte Solutions in the Water Supercooled Regime: Signatures of the Liquid-Liquid Water Phase Transition

TL;DR

This study uses molecular dynamics simulations to investigate how electrolyte solutions in supercooled water reveal anomalies linked to a hypothesized liquid-liquid phase transition, affecting thermal and thermoelectric properties.

Contribution

It demonstrates that electrolyte solutions preserve water's supercooled anomalies and links thermal transport behavior to structural water phase changes.

Findings

Thermal conductivity shows a minimum near 220 K at low salt concentrations.

Soret coefficient changes from thermophobic to thermophilic with cooling.

Seebeck coefficient reversals indicate structural transformations.

Abstract

Water exhibits remarkable anomalies when supercooled, attributed to a hypothesized liquid-liquid phase transition (LLPT) between low-density (LDL) and high-density (HDL) liquid phases. Using non-equilibrium molecular dynamics simulations, we explore thermal transport and coupled effects in supercooled NaCl and LiCl solutions (1-4 m, 200-300 K). At 1 m, thermal conductivity exhibits a pronounced minimum near 220 K, coinciding with maxima in isothermal compressibility and minima in the speed of sound, both of which are signatures of critical fluctuations. The anomalies progressively diminish with increasing salt concentration and vanish at 4 m, suggesting suppression of the LLPT. The Soret coefficient exhibits a striking behavior. Initially thermophobic at high temperatures (> 280 K), becoming thermophilic upon cooling, then reverting to thermophobic below 220 K. This behavior correlates with structural changes in the hydrogen-bond network of water. Specifically, we find that electrolyte solutions dominated by HDL structures, which are characterized by lower tetrahedral order, exhibit thermophobic behavior, whereas thermodynamic states dominated by LDL structures, with higher tetrahedral order, display thermophilic behavior. Furthermore, Seebeck coefficients exhibit sign reversals near 220-230 K, highlighting the thermoelectric sensitivity to structural transformations and temperature. These findings establish thermal transport as a sensitive probe of supercooled water, revealing that electrolyte solutions preserve the water's anomalies deep into the supercooled regime.

Tables

Table 1. Table 1: Average tetrahedral order parameter ξ ¯ \overline{\xi} , percentage of HDL area fraction (HDL), and Soret coefficient s T s_{T} for NaCl (top) and LiCl (bottom) aqueous solutions at 1, 2, and 4 m from 200–300 K temperature. ξ ¯ \overline{\xi} and HDL area fraction were obtained from tetrahedral order parameter distributions. The Soret coefficients are from the NEMD simulations with standard deviations from replicas reported in parentheses. Bold entries indicate the local minima of s T s_{T} and matched HDL% for the correlation analysis.

	1 m			2 m			4 m
T	$\overline{\xi }$	HDL	$s_T$	$\overline{\xi }$	HDL	$s_T$	$\overline{\xi }$	HDL	$s_T$
$\left[K\right]$		$[\%]$	$\left[{10}^{-3}{K}^{-1}\right]$		$[\%]$	$\left[{10}^{-3}{K}^{-1}\right]$		$[\%]$	$\left[{10}^{-3}{K}^{-1}\right]$
NaCl
200	0.773	16.5	4.187(0.8)	0.708	24.2	-0.040(0.8)	0.608	48.5	-3.549(0.6)
220	0.738	21.1	-0.767(0.6)	0.683	27.1	-4.018(0.4)	0.592	48.8	-2.325(0.5)
240	0.707	24.7	-3.129(0.4)	0.661	31.6	-1.726(0.5)	0.578	54.9	-0.855(0.3)
260	0.683	28.9	-0.747(0.4)	0.637	33.3	0.122(0.5)	0.564	68.1	0.193(0.3)
280	0.654	35.3	1.266(0.3)	0.615	40.3	1.221(0.2)	0.549	72.1	1.122(0.1)
300	0.628	34.3	2.316(0.3)	0.594	51.4	2.120(0.1)	0.534	75.5	1.850(0.1)
LiCl
200	0.766	15.7	0.449(0.5)	0.719	22.4	0.035(0.7)	0.625	36.9	-0.186(0.6)
220	0.750	18.3	-1.479(0.4)	0.700	23.8	-1.529(0.4)	0.606	40.9	-1.000(0.3)
240	0.717	22.9	-2.407(0.5)	0.672	27.2	-2.513(0.4)	0.586	45.5	-1.450(0.3)
260	0.688	27.6	-2.499(0.5)	0.646	33.9	-1.785(0.4)	0.565	49.4	-1.442(0.3)
280	0.658	32.8	-0.937(0.3)	0.620	35.1	-0.991(0.3)	0.545	57.4	-0.914(0.2)
300	0.632	32.8	0.199(0.2)	0.597	42.6	-0.042(0.2)	0.524	62.3	-0.297(0.2)

Table 2. Table 2: NEMD simulation results for NaCl and LiCl aqueous solutions at 1, 2, and 4 m. Columns list the mean temperature T T of the simulation cell, the mean temperatures of the hot and cold reservoirs T hot T_{\mathrm{hot}} and T cold T_{\mathrm{cold}} , density ρ \rho , pressure P P , thermal conductivity λ \lambda , and Seebeck coefficient S E S_{E} . Values are time averages over steady-state windows, and numbers in parentheses denote standard deviations across independent replicas.

T	T${}_{\text{hot}}$	T${}_{\text{cold}}$	$\rho$	$P$	$\lambda$	$S_E$
[K]	[K]	[K]	[g/cm3]	[atm]	[W m-1 K-1]	[mV K-1]
1m NaCl
200.2	249.9	150.0	1.009	8.9(22.2)	0.786(0.016)	2.233(0.627)
220.3	270.1	170.1	1.028	-4.2(9.7)	0.777(0.013)	0.294(0.555)
239.9	291.0	190.0	1.041	12.5(3.1)	0.792(0.010)	-0.299(0.458)
260.2	310.1	210.1	1.044	10.6(3.9)	0.808(0.009)	-0.386(0.276)
280.0	329.8	229.9	1.042	4.6(1.8)	0.826(0.009)	-0.422(0.131)
300.0	350.8	249.9	1.036	6.5(2.3)	0.840(0.004)	-0.403(0.056)
2m NaCl
200.1	250.1	150.0	1.066	10.2(25.5)	0.786(0.009)	0.853(0.573)
220.4	270.8	170.0	1.078	0.3(7.6)	0.787(0.008)	-0.065(0.527)
240.2	290.3	189.9	1.086	4.9(4.7)	0.791(0.006)	-0.270(0.468)
260.2	310.5	210.0	1.085	3.9(3.7)	0.803(0.007)	-0.470(0.202)
280.1	329.2	229.9	1.080	-1.5(1.2)	0.815(0.008)	-0.497(0.251)
300.0	350.3	249.9	1.067	8.5(1.8)	0.825(0.001)	-0.446(0.131)
4m NaCl
200.0	250.1	149.9	1.166	-88.6(34.0)	0.767(0.014)	0.741(0.665)
220.1	270.2	170.1	1.169	-72.0(7.1)	0.771(0.015)	0.025(0.456)
240.1	290.1	189.9	1.166	-63.2(4.9)	0.778(0.007)	-0.343(0.154)
260.1	310.3	209.7	1.158	-66.9(3.0)	0.786(0.007)	-0.414(0.230)
280.1	330.8	229.8	1.147	-66.2(1.0)	0.792(0.007)	-0.461(0.235)
299.8	350.3	249.9	1.134	-33.4(2.2)	0.801(0.004)	-0.491(0.133)
1m LiCl
200.3	250.5	149.9	0.985	38.7(16.1)	0.785(0.011)	1.606(0.733)
220.6	270.8	169.9	1.005	7.0(16.8)	0.775(0.010)	0.067(0.533)
240.3	290.8	190.1	1.015	-7.9(4.3)	0.783(0.008)	-0.437(0.281)
260.3	310.6	210.1	1.022	-2.3(2.3)	0.802(0.005)	-0.545(0.224)
280.2	330.9	229.9	1.023	1.1(1.5)	0.814(0.005)	-0.521(0.275)
300.1	350.2	249.9	1.018	1.8(0.8)	0.835(0.004)	-0.559(0.149)
2m LiCl
200.0	250.1	149.6	1.019	76.2(17.8)	0.780(0.009)	0.412(0.728)
220.5	270.6	169.9	1.036	-12.8(18.6)	0.780(0.008)	-0.245(0.482)
240.2	290.7	190.1	1.042	-9.4(2.7)	0.781(0.006)	-0.384(0.369)
260.2	311.0	210.0	1.046	-5.3(3.5)	0.794(0.006)	-0.300(0.262)
280.2	330.7	230.2	1.044	-1.7(1.0)	0.808(0.004)	-0.318(0.255)
300.1	350.5	249.9	1.039	4.3(1.1)	0.821(0.003)	-0.345(0.159)
4m LiCl
200.0	250.0	149.5	1.082	-49.8(37.6)	0.760(0.010)	0.269(0.777)
220.1	270.0	169.9	1.088	-7.3(14.4)	0.762(0.007)	-0.210(0.479)
240.2	290.6	190.2	1.089	-3.5(4.5)	0.769(0.006)	-0.256(0.257)
260.1	310.6	209.7	1.087	-0.2(2.1)	0.778(0.007)	-0.345(0.294)
280.1	330.9	229.9	1.084	0.4(0.9)	0.784(0.006)	-0.347(0.299)
300.0	350.5	250.0	1.078	-2.2(0.9)	0.796(0.003)	-0.464(0.176)

Table 3. Table 3: Thermophysical properties (density ρ \rho , isobaric thermal expansion coefficient α \alpha , isothermal compressibility κ T \kappa_{T} , adiabatic speed of sound c s c_{s} ) and structural properties (coordination numbers between the ions and oxygen in water, with cutoff radius r r at the first minimum of the ion–O g ​ ( r ) g(r) ) for NaCl and LiCl aqueous solutions at 1, 2, and 4 m over 200–300 K from equilibrium NPT simulations. Values are means and uncertainties in parentheses are standard deviations across independent replicas.

T	$\rho$	$\alpha$	${\kappa }_T$	$c_s$	n${}_{\text{Na–O}}$	$r_{\text{Na–O}}$	n${}_{\text{Cl–O}}$	$r_{\text{Cl–O}}$
[K]	[g/cm3]	[$\times {10}^{-4}$ K-1]	[$\times {10}^{-10}$ m3J-1]	[m/s]		[nm]		[nm]
1m NaCl
200	1.009	-2.54(0.58)	4.11(0.29)	1565.4(36.0)	5.68	0.31	5.76	0.36
220	1.028	-6.34(0.26)	5.79(0.25)	1301.5(31.4)	5.67	0.31	5.78	0.36
240	1.041	-3.96(0.14)	5.05(0.19)	1391.1(19.4)	5.63	0.31	5.79	0.36
260	1.044	0.79(0.18)	4.50(0.05)	1476.3(10.3)	5.61	0.32	5.80	0.36
280	1.042	2.15(0.13)	4.24(0.07)	1518.2(10.9)	5.57	0.32	5.80	0.36
300	1.036	3.71(0.16)	4.11(0.06)	1556.4(9.4)	5.50	0.31	5.80	0.36
2m NaCl
200	1.066	-8.04(0.56)	4.85(0.37)	1393.3(27.3)	5.78	0.32	5.82	0.36
210	1.070	-5.76(0.51)	5.40(0.27)	1335.6(27.6)	5.75	0.32	5.82	0.36
220	1.078	-5.29(0.32)	4.83(0.22)	1386.9(17.1)	5.73	0.31	5.82	0.36
240	1.086	-1.50(0.20)	4.17(0.16)	1486.1(9.9)	5.68	0.32	5.82	0.36
260	1.085	1.83(0.16)	3.95(0.14)	1532.8(10.5)	5.62	0.32	5.89	0.36
280	1.080	3.28(0.07)	3.72(0.04)	1581.6(9.2)	5.56	0.32	5.90	0.36
300	1.067	4.27(0.07)	3.72(0.04)	1614.5(9.8)	5.51	0.32	5.93	0.37
4m NaCl
200	1.166	-2.87(0.65)	4.16(0.26)	1437.5(27.9)	5.85	0.32	5.87	0.36
220	1.169	-0.64(0.32)	3.85(0.20)	1515.4(29.3)	5.78	0.32	5.89	0.36
240	1.166	1.67(0.27)	3.40(0.13)	1607.7(20.5)	5.70	0.32	5.90	0.36
260	1.158	2.95(0.19)	3.18(0.03)	1663.9(12.8)	5.61	0.32	5.91	0.36
280	1.147	4.01(0.09)	3.15(0.03)	1689.7(10.7)	5.51	0.31	5.93	0.36
300	1.134	4.77(0.09)	3.14(0.01)	1707.8(7.1)	5.47	0.32	5.94	0.36
T	$\rho$	$\alpha$	${\kappa }_T$	$c_s$	n${}_{\mathrm{Li}\text{–}O}$	$r_{\mathrm{Li}\text{–}O}$	n${}_{\text{Cl–O}}$	$r_{\text{Cl–O}}$
[K]	[g/cm3]	[$\times {10}^{-4}$ K-1]	[$\times {10}^{-10}$ m3J-1]	[m/s]		[nm]		[nm]
1m LiCl
200	0.985	-6.33(0.44)	4.51(0.42)	1503.0(42.4)	4.00	0.24	5.62	0.36
220	1.005	-7.32(0.25)	5.02(0.20)	1411.0(38.9)	4.00	0.24	5.71	0.36
240	1.015	-3.89(0.15)	5.01(0.15)	1428.6(30.1)	4.00	0.25	5.71	0.36
260	1.022	-1.72(0.15)	4.43(0.07)	1493.2(14.9)	4.00	0.26	5.76	0.36
280	1.023	1.56(0.29)	4.31(0.05)	1524.4(17.9)	4.00	0.26	5.82	0.37
300	1.018	3.50(0.19)	4.19(0.05)	1540.5(12.5)	4.00	0.26	5.85	0.37
2m LiCl
200	1.019	-6.50(0.37)	4.01(0.20)	1563.4(31.3)	4.00	0.23	5.64	0.36
220	1.036	-5.19(0.45)	4.19(0.19)	1549.7(20.1)	4.00	0.24	5.75	0.36
240	1.042	-3.04(0.18)	4.02(0.11)	1551.7(17.1)	4.00	0.25	5.75	0.36
260	1.046	-0.68(0.18)	3.93(0.07)	1569.7(12.4)	4.00	0.25	5.76	0.36
280	1.044	2.19(0.08)	3.84(0.07)	1584.8(12.0)	4.00	0.26	5.83	0.37
300	1.039	3.79(0.12)	3.78(0.06)	1618.3(12.1)	4.00	0.28	5.90	0.37
4m LiCl
200	1.082	-2.63(0.51)	3.61(0.18)	1597.6(31.5)	4.00	0.24	5.74	0.36
220	1.088	-0.75(0.43)	3.48(0.14)	1634.4(24.7)	4.00	0.24	5.73	0.36
240	1.089	-0.24(0.37)	3.26(0.09)	1678.7(15.8)	4.00	0.25	5.74	0.36
260	1.087	0.75(0.12)	3.24(0.06)	1694.6(12.0)	4.00	0.25	5.79	0.36
280	1.084	2.55(0.14)	3.18(0.05)	1715.3(11.5)	4.00	0.26	5.81	0.36
300	1.078	4.03(0.11)	3.14(0.06)	1728.4(11.6)	4.00	0.26	5.89	0.37

Equations

$$\lambda (z)=-\frac{J_q}{\nabla T(z)}$$

$$J_q=\frac{\dot{Q}}{2A},$$

$$s_T(T)=-\frac{1}{x_1(T)(1-x_1(T))}{\left(\frac{\nabla x_1(T)}{\nabla T(T)}\right)}_{J_1=0}=-\frac{1}{b(T)}{\left(\frac{db(T)}{dT}\right)}_{J_1=0},$$

$$\xi =1-\frac{3}{8}\sum _{j=1}^3\sum _{k=j+1}^4{\left(\cos {\psi }_{jk}+\frac{1}{3}\right)}^2$$

$$\rho_{q}(z)=\frac{1}{A}\Biggl{\langle}\sum_{i=1}^{N_{charges}}q_{i}\delta(z-z_{i})\Biggr{\rangle},$$

$$E_z(z)=\frac{1}{{\epsilon }_0}{\int }_0^z{\rho }_q(z^{\prime })d{z}^{\prime },$$

$$S_E(z)=\frac{E_z(z)}{\nabla T(z)},$$

$$c_s=\sqrt{\frac{C_P}{C_V}\frac{1}{\rho {\kappa }_T}}.$$

$$\varphi (z)=-{\int }_0^z{E}_z(z^{\prime })d{z}^{\prime }.$$

Full text

Keywords: Supercooled water, Liquid-liquid phase transition, Thermal conductivity, Soret effect Molecular dynamics simulation, Electrolyte solutions

I Introduction

Water, being the most abundant liquid on Earth, plays a central role in various natural and technological processes within the physical, chemical, and biological sciences.Franks (2000) Unlike most liquids, water exhibits a density maximum and displays anomalous behaviors, such as increased isothermal compressibility and heat capacity when cooled.Eisenberg and Kauzmann (2005); Debenedetti (1996) These anomalies become even more pronounced in the supercooled regime,Speedy and Angell (1976); Angell, Sichina, and Oguni (1982); Debenedetti and Stanley (2003) raising fundamental questions about the underlying mechanisms governing water’s structure and dynamics.

Several hypotheses have been proposed to explain these anomalies, with one of the most compelling being the liquid-liquid phase transition (LLPT) hypothesis. This hypothesis suggests that there is a transition between a low-density liquid (LDL) phase, which exhibits strong tetrahedral order, and a high-density liquid (HDL) phase, characterized by weaker local order.Poole et al. (1992) This transition is believed to terminate at a liquid–liquid critical point. The LLPT hypothesis attributes the anomalous thermodynamic behavior of supercooled water to critical fluctuations occurring as the liquid approaches the second critical point.Gallo et al. (2016)

Recent simulation studies employing empirical and ab initio trained machine learning-based models have provided increasing evidence supporting the LLPT hypothesis.Palmer et al. (2014); Debenedetti, Sciortino, and Zerze (2020); Donkor et al. (2024); Malek et al. (2025); Sciortino et al. (2025) Furthermore, experimental investigations, although technically challenging due to the rapid crystallization of water, have also yielded tentative support for the existence of a second critical point.Mishima and Stanley (1998); Gallo et al. (2016); Kim et al. (2017, 2020); Amann-Winkel et al. (2023)

The studies mentioned earlier clearly demonstrate the existence of thermodynamic anomalies in pure water. However, it is crucial to understand how these anomalies behave in aqueous solutions, both for fundamental research and practical applications. This understanding is important because, in both natural and technological contexts, water rarely exists on its own and usually acts as a solvent. In many cases, water can be more readily supercooled in solution than in its pure form.Kanno and Angell (1977) Simulation studies of NaCl and LiCl aqueous solutions using various force fields support the existence of a liquid–liquid critical point (LLCP), which shifts to higher temperatures and lower pressures upon salt addition.Corradini, Rovere, and Gallo (2010); Corradini and Gallo (2011); Perin and Gallo (2023) Notably, the LLCP has been reported to disappear in concentrated LiCl solutions (2.0 m).Perin and Gallo (2023) These findings motivate further investigation into how ion type and concentration influence the LLCP location, and whether these structural changes are connected to the anomalies in supercooled aqueous electrolyte solutions.

Previous investigations into the thermal transport properties of supercooled pure water have reported a minimum in thermal conductivity Biddle et al. (2013); Kumar and Stanley (2011); Bresme et al. (2014); Zhao and Bresme (2024). This minimum has been connected to various thermodynamic anomalies, such as maximum isothermal compressibility and minimum adiabatic speed of sound. These phenomena are believed to be associated with the hypothesized liquid-liquid critical point (LLCP).

The thermal conductivity quantifies the magnitude of the heat flux emerging from an imposed temperature gradient. Thermal gradients in pure polar fluids induce the fluid polarization, an effect known as the thermal polarization (TP) effectBresme et al. (2008); Chapman and Bresme (2022). In electrolyte solutions, thermal gradients give rise to thermodiffusion (Ludwig-Soret effect)Ludwig (1856); Soret (1879) and thermoelectricity (Seebeck effect)Seebeck (1826). These non-equilibrium effectsde Groot and Mazur. (1984) involve coupling between mass and heat fluxes, resulting in the development of concentration gradients (Soret) in fluid mixtures or electrostatic potential differences (Seebeck) in metals, as well as in both aqueous and non-aqueous solutions.Bonetti et al. (2015, 2011); Di Lecce and Bresme (2018); Nickel et al. (2024) Indeed the Seebeck effect has gained increasing attention for enabling electric and electrochemical devices to harvest low-grade heat through temperature-induced voltages, with ionic systems offering promising energy conversion.

The Soret effect is an intriguing transport phenomenon,Wiegand (2004); Köhler and Morozov (2016) that is increasingly being used in practical applications as a powerful tool for assessing protein-ligand binding affinitiesWienken et al. (2010); Seidel et al. (2013); Niether and Wiegand (2019). However, we still lack a comprehensive theory to predict whether a solute will migrate towards hot (thermophilic) or cold (thermophobic) regions. Additionally, the Soret coefficient often exhibits an exotic behavior that poses challenges to microscopic interpretation. For instance, minima in the Soret coefficient have been observed in both experiments and simulations of electrolyte solutions, revealing a strong correlation between these minima and the overall thermodiffusion response.Colombani, Bert, and Dupuy-Philon (1999); Di Lecce, Albrecht, and Bresme (2017a); Rudani, Briels, and Wiegand (2025) Recent simulation studies have also demonstrated a relationship between thermophobicity in electrolyte solutions and the disruption of the tetrahedral structure of water.Zhao and Bresme (2025) This finding aligns with earlier conclusions drawn from experimental analyses of aqueous solutions.Niether and Wiegand (2019)

In this work, we employ non-equilibrium molecular dynamics (NEMD) simulations to explore the thermal transport properties of NaCl and LiCl aqueous solutions within the supercooled water regime. We investigate the thermophysical and structural properties in these supercooled states, with a particular focus on isothermal compressibility due to its relevance to the hypothesized liquid-liquid phase transition (LLPT).

We demonstrate that the thermodiffusive response of aqueous alkali halides is quite complex, featuring multiple thermophobic and thermophilic transitions. This complexity is evident from the presence of minima in the Soret coefficient and a prevalent thermophobic response observed at both high and very low temperatures. In contrast, a thermophilic response occurs at intermediate temperatures. Furthermore, we identify a minimum in thermal conductivity at a concentration of 1 m, which disappears in highly concentrated solutions, 4 m. We also discuss the potential thermodynamic and structural origins of these observations.

II Methods

Water molecules were modeled using the TIP4P/2005 model,Abascal and Vega (2005) a non-polarizable, rigid water model widely recognized for accurately reproducing the properties of supercooled water and various ice polymorphs.Vega and Abascal (2011); Vega et al. (2008) Ions were simulated using the Madrid 2019 force field,Zeron, Abascal, and Vega (2019), which is specifically parameterized for compatibility with TIP4P/2005 and employs scaled ionic charges to capture the thermodynamic and structural properties of electrolyte solutions. This force field has been shown to reliably reproduce experimental thermal transport behavior of NaCl and LiCl under ambient conditionsBresme and Vasey (2024), and it has also been applied in simulations of supercooled electrolyte solutions.Perin and Gallo (2023); Sedano, Vega, and Noya (2024)

We used boundary-driven non-equilibrium molecular dynamics (NEMD) simulations, performed with LAMMPS (version from 29 Aug 2024),Plimpton (1995) to investigate the thermal transport and the induced coupling effects of the alkali halide aqueous solutions. We applied thermostats to designated regions at the center and edges of the box, maintaining predefined hot ($T_{\text{hot}}$) and cold ($T_{\text{cold}}$) temperatures (see Figure 1). The thermostatting regions were 8 Å thick, and the simulation box was configured with dimensions in the ratio $\{L_{x}:L_{y}:L_{z}\}=\{1:1:3\}$, where $L_{x}$ was approximately $40\,~\text{\AA }$, depending on the density.

Langevin thermostatsSchneider and Stoll (1978); Brünger, Brooks, and Karplus (1984) were employed to control the temperature of molecules and ions at each timestep (1 fs) with a coupling constant of 500 fs. We also removed the center of mass momentum of the entire system at every timestep. Electrostatic interactions were computed using the Particle-Particle Particle-Mesh (PPPM) method, Hockney and Eastwood (1989) with a relative error of 10*-5*. The PPPM method was likewise employed to handle long-range dispersion interactions,Isele-Holder et al. (2013); Isele-Holder, Mitchell, and Ismail (2012) with accuracy thresholds set to 10*-4* kcal/(mol Å) for the real space component and 2$\times$10*-4* kcal/(mol Å) for the k-space component.

A representative simulation box contained 6663 water molecules modeled using the TIP4P/2005 potential,Abascal and Vega (2005) along with 120, 240, or 480 ion pairs, corresponding to concentrations of 1, 2, and 4 m, respectively. Ion-ion and ion-water interactions were described using the Madrid-2019 force field.Zeron, Abascal, and Vega (2019) Prior to the NEMD runs, the system was pre-equilibrated for 5 ns in the isothermal-isobaric (NPT) ensemble at a pressure of 1 bar and a temperature set to the average $T_{\text{hot}}$ and $T_{\text{cold}}$. The equilibration process employed a Nosé-Hoover thermostat and barostat with relaxation times of 100 fs and 1000 fs, respectively.Nosé (1984); Hoover (1985) Then the thermostats were activated, and the NEMD simulations were conducted under constant volume conditions.

We carried out long simulations for supercooled systems. The simulations lasted at least 200 ns, discarding the initial 150 ns to ensure the system had reached a steady state (see Figure 8 in Appendix A). Following this initialization period, we computed the local temperature and concentration profiles by averaging over 200 spatial bins along the $z$-axis, corresponding to a bin width of approximately $0.6\,~\text{\AA }$. These profiles were further averaged over a production run of at least 50 ns across 10 independent replicas, each initialized with a different random seed for the thermostats, yielding a total sampling time of 0.5 $\mu$s. The trajectory from at least the final 10 ns of the production run was saved every 1000 steps, resulting in a total of $10^{5}$ configurations used to compute tetrahedral order parameters and electrostatic potential profiles.

The thermal conductivity, $\lambda(z)$, along the direction of the local thermal gradient $\nabla T(z)$, was computed with Fourier’s Law,

[TABLE]

Here, $J_{q}$ represents the heat flux, which was determined using the continuity equation:

[TABLE]

where $\dot{Q}$ denotes the rate of energy exchanged at the hot and cold thermostats, $A$ represents the cross-sectional area of the simulation box perpendicular to the direction of heat flow, and the factor of 2 accounts for the presence of the two heat fluxes within the simulation cell. Our approach ensures excellent energy conservation, with equal amounts of energy transferred at the cold and hot thermostats within computational accuracy (see Figure 9 in Appendix B).

The Soret coefficient $s_{T}$ was determined from the local molar fraction profiles obtained in NEMD production runs:

[TABLE]

where $x_{1}(T)$ denotes the molar fraction of salt at the local temperature $T$, respectively, and $b(T)$ is the corresponding solute molality. Typical local concentration profiles used for the Soret coefficient calculations are shown in Appendix C (see Figure 10).

We also calculated the tetrahedral order parameter $\xi$ to investigate structural changes in the supercooled electrolyte solutions, using the following expression:Errington and Debenedetti (2001)

[TABLE]

where $\psi_{jk}$ is the angle between the lines connecting a central water molecule and its four nearest oxygen neighbors.

To investigate the thermoelectric effect in electrolyte solutions, we calculated the charge density $\rho_{q}(z)$,

[TABLE]

where the angle brackets denote a time-averaged binning of the sum over all charged sites. In this expression, $q_{i}$ is the charge of site $i$, and $z_{i}$ is its coordinate along the $z$-axis. In the system studied, the charged sites comprise the hydrogen ($q_{H}=0.5564e$) and dummy atom ($q_{M}=-1.1128e$) of the TIP4P/2005 water model.Abascal and Vega (2005) Additional contributions come from the scaled charges of the ions, modeled using the Madrid-2019 potential: ($q_{Na^{+}}=q_{Li^{+}}=0.85e$) and ($q_{Cl^{-}}=-0.85e$).Zeron, Abascal, and Vega (2019) The electric field was obtained by integrating the charge density:

[TABLE]

where $\varepsilon_{0}$ is the vacuum permittivity. To quantify the strength of the thermoelectricity response, we calculated the Seebeck coefficient $S_{E}$,

[TABLE]

which represents the ratio of the $z$ component of the electrostatic field to the thermal gradient. The detailed procedure for calculating the Seebeck coefficient from the local electrostatic potential profiles is provided in Appendix D.

To investigate the thermophysical properties of the electrolyte solutions, we performed equilibrium molecular dynamics simulations in the isothermal–isobaric (NPT) ensemble using GROMACS (vs. Jan 2024).Bekker et al. (1993); Abraham et al. (2024) Each system consisted of a cubic simulation box containing 926 water molecules and 17, 34, or 67 ion pairs, corresponding to 1.0, 2.0, and 4.0 m concentrations, respectively. The production runs were carried out for at least 200 ns with a timestep of 1 fs to ensure adequate sampling. The fluctuation-derived thermodynamic properties, including the isobaric thermal expansion coefficient ($\alpha$), isothermal compressibility ($\kappa_{T}$), the isobaric ($C_{P}$) and isochoric ($C_{V}$) heat capacities, were computed using the following fluctuation relations:Allen and Tildesley (1987) $\langle\delta V^{2}\rangle_{NPT}=Vk_{B}T\kappa_{T},\;C_{P}=\left<\delta H^{2}\right>/k_{B}T^{2},\;C_{V}=C_{P}-VT\alpha^{2}/\kappa_{T}.$ In these expressions, $V$ is the system volume, $H$ is the enthalpy, and $k_{B}$ is the Boltzmann constant. The fluctuation $\delta\mathcal{P}=\mathcal{P}-\mathcal{P}_{ens}$ denotes the deviation of a thermodynamic property $\mathcal{P}$ from its ensemble average $\mathcal{P}_{\text{ens}}$ in the NPT ensemble. The isentropic speed of sound was then calculated from these fluctuation properties along with the mass density $\rho$:

[TABLE]

III Results

Figure 2 illustrates the equations of state for NaCl and LiCl aqueous solutions at concentrations of 1, 2, and 4 molal (m), calculated along the standard pressure isobar (1 atm) using local temperatures and densities obtained from Non-Equilibrium Molecular Dynamics (NEMD) simulations. The results demonstrate a strong agreement with our own equilibrium NPT simulation data, as well as with previous molecular dynamics simulations Sedano et al. (2022). This consistency between our equilibrium and non-equilibrium results, along with earlier studies, reinforces the validity of the simulation conditions applied in our NEMD computations. Our results indicate that the systems achieve local equilibrium under thermal gradients of approximately 1 K/Å. In line with findings reported in earlier studies, our results show that the addition of salt causes a shift in the temperature of the density maximum (TMD) to lower temperatures. As discussed by Sedano et al. (2022), the predictions made using the TIP4P/2005 - Madrid-2019 force field align closely with the experimental TMDs of LiCl and NaCl within the 0-2 m concentration range.

Figure 3 shows the thermal conductivities of NaCl and LiCl aqueous solutions as a function of temperature within the range of 200–300 K, at concentrations varying from 1 to 4 m. The thermal conductivity of these solutions decreases with temperature. However, for salt concentration 1 m, a minimum in the thermal conductivity is observed at lower temperatures, around 220 K. The results also show that the thermal conductivity decreases as the density of the fluid increases in a wide temperature range, indicating an anomalous behavior in the thermal response of the solution. Advancing the discussion below, this reduction is associated with a decrease in the speed of sound in the solution throughout most of the supercooled region examined here.

The thermal conductivity of NaCl and LiCl solutions is very similar across the studied temperature and concentration ranges. This behavior is consistent with previous simulation studies, which reported that the thermal conductivity of alkali halide solutions is more strongly influenced by anion mass and composition than by cation type. Bresme and Vasey (2024) At 300 K, both 1 m NaCl and LiCl solutions exhibit conductivities close to that of pure water, while at higher concentrations a pronounced reduction is observed.

In both NaCl and LiCl solutions at 1 m concentration, the thermal conductivity exhibits a local minimum near 220 K, similar to what is observed in pure TIP4P/2005 water.Zhao and Bresme (2024) However, this minimum disappears at higher concentrations (2 m and 4 m) indicating that salt addition suppresses the thermal conductivity anomaly.

The minimum in thermal conductivity suggests an underlying thermodynamic origin. To understand this connection, we examine the isothermal compressibility, which is expected to influence thermal transport due to its relationship to the speed of sound. The compressibility was computed using equilibrium simulations conducted in the NPT ensemble and was used to calculate the isentropic (adiabatic) speed of sound. Our findings on isothermal compressibility are in good agreement with both experimental data and previous simulation results, under both ambient and supercooled conditions Rogers and Pitzer (1982); Sedano, Vega, and Noya (2024); Bresme and Vasey (2024) (see Figure 4).

The isothermal compressibility of the NaCl and LiCl 1 m and 2 m solutions exhibits prominent maxima, resembling the behavior observed for pure water in earlier studies Gallo et al. (2016). However, the general trend indicates that these maxima shift to lower temperatures and decrease in magnitude. This behavior is consistent with the findings reported in previous works Perin and Gallo (2023); Sedano, Vega, and Noya (2024), and it has significant implications for the behavior of the speed of sound (see Figure 4). The calculated speed of sound for NaCl aligns well with available experimental data Kleis and Sanchez (1990). It decreases with temperature in the moderate undercooling regime (T $>$ 240 K) and in our simulations it features minima for the NaCl and LiCl 1 m solutions at temperatures corresponding to the points where compressibility exhibits a maximum. These minima coincide with the region where a minimum in thermal conductivity is also observed (see Figure 3). We note that the influence of the speed of sound on thermal conductivity was previously noted by Bridgman in his empirical equation Bridgman (1923).

At higher concentrations, specifically 4 m, the speed of sound does not feature any distinct minima. The lack of clear minima within the temperature range studied may be attributed to a shift in the temperature of these extrema to values lower than those included in our investigation, as well as a further attenuation of the anomalies in more concentrated solutions. These findings collectively support the view that thermal transport in supercooled alkali halide solutions remains influenced by underlying thermodynamic anomalies. This interpretation suggests that the minimum in thermal conductivity is aligned with the presence of a liquid-liquid critical point (LLCP), which is evident in the 1 m solution. However, at higher concentrations, these anomalies are significantly diminished or no longer detectable, consistent with previous research indicating that the LLCP is suppressed or completely disappears in concentrated lithium chloride (LiCl) solutions.Perin and Gallo (2023)

Our analysis of the speed of sound indicates a correlation with the lack of minima in the thermal conductivity of highly concentrated solutions. This suggests that the reduction in thermal conductivity observed across all salts and concentrations, along with the minima found in 1 m solutions, is related to the temperature dependence of the speed of sound. Ultimately, this behavior arises from the temperature dependence of the isothermal compressibility of the solutions in the supercooled regime.

Having established the thermodynamic signatures of structural transitions, we now examine thermodiffusive transport, which provides insight into the relationship between thermal gradients and molecular mobility. In Figure 5, we present the Soret coefficients as a function of temperature for aqueous NaCl and LiCl solutions at various concentrations. At 300 K, our results align with previous simulation studies, accurately reflecting both the magnitude and temperature dependence of the Soret effect.Bresme and Vasey (2024) Notably, we observe a transition in the Soret coefficient from thermophobic behavior (at high temperatures) to thermophilic behavior (at low temperatures). This transition occurs for LiCl at approximately 300 K and for NaCl around 260 K, regardless of the concentration studied. This sign reversal has also been documented in experimental work and simulations using different force fields.Gaeta et al. (1982); Römer et al. (2013); Lee et al. (2024); Di Lecce, Albrecht, and Bresme (2017b); Di Lecce and Bresme (2019)

A notable feature is the crossing of the Soret coefficient curves for different concentrations near 280 K. Above this temperature, the Soret coefficient decreases with increasing concentration, whereas below it, the trend is reversed. A similar crossing behavior was observed in experiments at approximately 315 K,Römer et al. (2013) as well as in simulations of LiCl solutions using the SPC/E water model, where the crossing occurred around 260–280 K.Di Lecce, Albrecht, and Bresme (2017b, a) This temperature interval is close to the value obtained here. These previous investigations attributed the concentration dependence to competing mechanisms: at high temperatures, entropy predominates, while at low temperatures, stronger hydrogen bonding enhances ion–water interactions, leading to an accumulation of salt in the heated region. Our results replicate these qualitative trends and extend them into the supercooled regime.

In the deeply supercooled regime, for both NaCl and LiCl (1 and 2 m), the Soret coefficient exhibits a distinct local minimum between 220 and 250 K (see Figure 5). This minimum shifts to lower temperatures as the concentration increases from 1 to 2 m. At a concentration of 4 m, NaCl does not exhibit any evidence of a minimum within the studied temperature range, while for LiCl, the minimum is not so clearly defined as compared with lower concentrations.

When comparing NaCl to LiCl, the Soret coefficient for NaCl displays a more significant change over a broader temperature range upon cooling, developing into a sharper minimum. Additionally, at temperatures close to 200 K, a second sign reversal occurs in the Soret coefficients of the 1 and 2 m solutions of both NaCl and LiCl. This transition is especially pronounced for the 1 m NaCl solution, where the Soret coefficient rises above $4\times 10^{-3}\,\mathrm{K^{-1}}$ at 200 K, indicating a strong thermophobic response. Generally, our results suggest that supercooled solutions become less thermophilic with increasing concentration in the temperature range of 230-280 K and less thermophobic for temperatures above 280 K. This indicates an overall weakening of the thermodiffusion response as salt concentration increases.

The complex behaviour of the Soret coefficient suggests changes in the hydrogen-bond network of water. To explore the structural origins of the observed thermodiffusion behavior, we examine the tetrahedral order parameter. Figure 6 presents the distribution of the tetrahedral order parameter. The emergence of a bimodal distribution is consistent with previous findings in supercooled pure water, reflecting the coexistence of distinct local environments.Errington and Debenedetti (2001); Kuo et al. (2021) In particular, the two dominant peaks correspond to more tetrahedrally ordered configurations associated with a low-density liquid–like (LDL) structures $(\xi\approx 0.8{-}1)$ and less ordered, high-density liquid-like (HDL) arrangements $(\xi\approx 0.4{-}0.6)$. A comparison across temperatures and concentrations indicates a tendency toward stronger tetrahedral ordering at lower temperatures and lower salt concentrations.

Interestingly, for all the LiCl solutions, an additional peak emerges at $\xi<0.2$, indicating the existence of highly disordered local environments. This multimodal distribution is absent in the case of NaCl, suggesting the emergence of heterogeneous structures in LiCl, which lead to the distortion of the water tetrahedral structure.

Due to its small radius, lithium produces a strong electrostatic field and is solvated by a stable hydration shell consisting of four to five water molecules, which involves a high heat of hydration Mähler and Persson (2012); Mason et al. (2015). Both classical simulations and ab initio computations predict a solvation shell made up of four water molecules arranged in a tetrahedral motif.Zeron, Abascal, and Vega (2019); Rempe et al. (2000) Experiments also indicate that the hydration number is relatively independent of ion concentration.Mason et al. (2015) These findings align with our calculations of the lithium-oxygen (Li-O) coordination number (CN) (see Table LABEL:table:NPT-data in Appendix G). Our simulations demonstrate that the CN is fairly consistent across different temperatures in the supercooled regime, which corroborates recent investigations Perin and Gallo (2023). In contrast, the coordination numbers for chloride-oxygen (Cl-O) and sodium-oxygen (Na-O) exhibit greater variability with respect to both temperature and concentration. Therefore, the observed increase in the height of the maximum across the 1-4 molality range (see Figure 6) is likely reflecting the larger presence of [Li(H2O)4]+ structures in the solution as the concentration increases.

We have analyzed the probability distributions in Figure 6 to quantify the relative fractions of HDL structures and their correlation with the Soret coefficients. The numerical data are listed in Table 1, and the detailed procedure for this calculation is provided in Appendix E. Between 240 and 300 K, the Soret coefficient for both NaCl and LiCl at all studied concentrations becomes more negative, i.e the solution becomes more thermophilic, with the ions migrating preferentially towards the hot region. This observation is consistent with previous discussions that greater water ordering enhances the thermophobic character of aqueous solutions.Niether and Wiegand (2019); Zhao and Bresme (2025) At around 240 K, where the 1 m and 2 m NaCl and LiCl solutions exhibit local minima in the Soret coefficient, the calculated area fraction of the HDL structure is $\sim$ 27%, indicating that thermophilicity is observed in general for larger populations of LDL the structure, where the coordination of water molecules with other water molecules features a stronger orientational structure and lower coordination. Indeed, we find that over 70% of the systems investigated in this work have an HDL area fraction below 30%. Similarly, we find that all the systems with a high content of HDL, $>65$ % have positive Soret coefficients, and are thermophobic.

As we approach the supercooled regime, we observe a rapid increase in the Soret coefficient. It becomes strongly positive, $\sim+4\times 10^{-3}$ K*-1*, for the 1 M NaCl solution and slightly positive, $\sim+0.5\times 10^{-3}$ K*-1* for the LiCl solution. These results suggest the presence of an additional competing mechanism influencing thermodiffusion in the deeply supercooled regime. It is possible that the order parameter we are examining is not sufficiently sensitive to capture other structural changes occurring in the supercooled aqueous solutions.

The relationship between water structure and thermodiffusion also applies to other transport phenomena. The Seebeck serves as an additional probe to examine how structural changes affect the thermal transport in the supercooled regime. (see Figure 7). Above 240 K, no significant differences are observed between NaCl and LiCl at any of the studied concentrations. Our results for 1 m NaCl at 300 K closely agree with a previous simulation study and fall within the reported error bars, despite the use of different force fields.Nickel et al. (2024) In the deeply supercooled region, we observe an inversion in the sign of the Seebeck coefficient around 220 K. Overall, our results indicate a significant Seebeck effect in the supercooled regime, with an enhancement of this effect in the deeply supercooled regime.

IV Conclusions

We conducted nonequilibrium molecular dynamics simulations to investigate thermal transport properties of NaCl and LiCl solutions in the supercooled water regime using TIP4P/2005 and Madrid-2019 force fields. Our results provide evidence supporting the liquid-liquid phase transition (LLPT) hypothesis in electrolyte solutions at moderate salt concentrations, $\sim$ 1 m.

At 1 m salt concentration, the thermal conductivity shows pronounced minima near 220 K. This temperature coincides with maxima in isothermal compressibility and minima in the speed of sound. These anomalies progressively weaken with increasing concentration and disappear at 4 m, consistent with salt-induced suppression of the hypothesized LLPTPerin and Gallo (2023). The correlation between thermal transport anomalies and thermodynamic signatures suggests these phenomena share a common structural origin.

The Soret effect reveals complex thermodiffusion behavior with multiple sign reversals. Solutions transition from thermophobic at high temperatures ($>$280 K) to thermophilic upon cooling, followed by a return to thermophobic behavior in the deeply supercooled regime ($<$220 K). Local minima in the Soret coefficient occur between 220-250 K, with NaCl solutions displaying sharper anomalies than LiCl. These transitions correlate with structural changes in the hydrogen bonding network, as evidenced by tetrahedral order parameter analysis, which indicates an interconversion between the HDL and LDL structures as temperature changes. Notably, LiCl solutions exhibit an additional population of highly disordered local environments ($\xi<$ 0.2), absent in NaCl, indicating that lithium ions induce stronger structural heterogeneity in the hydrogen-bond network, which may contribute to the weaker and less sharply defined thermodiffusion anomalies observed for LiCl compared to NaCl.

The Seebeck coefficients are on the order of 1 mV/K, and show sign inversions around 220-230 K, demonstrating thermoelectric sensitivity to structural transitions and confirming coupling between thermal and electrical responses.

This study establishes thermal transport measurements as sensitive probes of structural transitions in supercooled electrolyte solutions, offering new microscopic insights into water’s anomalous behavior. In particular, we demonstrate a correlation between the HDL/LDL structures and thermodiffusion. At high temperatures, when the HDL (lower tetrahedral order) structure is dominant, the solutions exhibit thermophobic responses. In contrast, they become thermophilic at temperatures where LDL, and therefore a stronger tetrahedral order, dominates.

Additional molecular simulations with different force fields and experimental validation at low temperatures and high salt concentrations would provide a route to test these theoretical predictions and advance our understanding of water’s phase behavior under extreme conditions.

DATA AVAILABILITY

The data that supports the findings of this study are available within the article.

Acknowledgements.

We acknowledge the ICL RCS High-Performance Computing facility and the UK Materials and Molecular Modelling Hub for computational resources, partially funded by the EPSRC (Grant Nos. EP/P020194/1 and EP/T022213/1.

AUTHOR DECLARATIONS

Conflict of Interest: The authors have no conflicts to disclose.

Appendix A Time Evolution and Equilibration at Low Temperatures

Figure 8 shows the time evolution of pressure and thermal conductivity at the lowest simulated temperature, 200 K. Similar to previous findings for supercooled pure water,Zhao and Bresme (2024) the alkali halide aqueous solutions display slow dynamics at low temperatures, requiring an extended period to reach equilibrium. Consequently, the initial 150 ns of the production run were excluded from the analysis.

Appendix B Energy Conservation Check

Figure 9 presents the energy exchanged between the cold and hot thermostatting regions during a representative 5 ns NEMD simulation. The consistency between the energy exchanged at the hot and cold thermostats confirms that the simulations feature good energy conservation.

Appendix C Representative Local Concentration Profiles and Soret Coefficient

Figure 10 shows typical plots of local molality concentration as a function of temperature from NEMD simulations of a 1 M NaCl solution. The Soret coefficient can be obtained from the temperature derivative of these concentration profiles according to Equation 3. The left and right panels correspond to simulations performed at the average temperatures 220 and 280 K, respectively. The concentration features a local minimum and maximum in salt concentration around 215 and 270 K. These extrema mark the approximate temperatures where the Soret coefficient changes sign, consistent with the trend observed in Figure 5 in the main text. Note that due to significant fluctuations, abrupt concentration jumps, and large uncertainties near the glass transition, particularly below 200 K, data in this low-temperature regime were excluded from the Soret coefficient calculations.

Appendix D Calculation of the Local Electrostatic Potential and Seebeck Coefficient

Figure 11 shows the dependence of electrostatic potential on temperature. The local electrostatic potential $\phi(z)$ was determined by first integrating the charge density to calculate the electric field, as described in Equation 6. The $z$-component of the electric field was then integrated to yield the electrostatic potential:

[TABLE]

The Seebeck coefficient was evaluated using the relation $S_{E}(z)=-\partial\phi(z)/\partial T(z)$. We fitted the electrostatic potential as a function of temperature with a 4th-order polynomial and evaluated the negative first derivative of the fit. Both panels in Figure 11 were calculated by analyzing trajectories spanning 5 ns, with the electrostatic potential value set to zero at the starting point of the integration.

At lower temperatures, as shown in the left panel, the electrostatic potential profiles display increased fluctuations compared to higher-temperature runs. This behavior is primarily attributed to the slow dynamics and reduced sampling efficiency in the low-temperature regime. To improve statistical reliability, the trajectory analysis was based on simulations extended to 25 ns with additional replicas. Nonetheless, the Seebeck coefficient at temperatures below 220 K still exhibits larger error bars, as shown in Figure 7. Despite this, a local minimum is observed around 220–230 K, indicating a sign reversal in the Seebeck effect near this temperature.

Appendix E Tetrahedral Order Parameter Distributions and Structural Fractionation

Figure 12 presents normalized representative probability densities $P(\xi)$ of the tetrahedral order parameter $\xi$, computed from NEMD trajectories using Equation 4. To quantify the relative populations of local structural motifs, the distributions are partitioned at the local minima between peaks, separating the less tetrahedrally coordinated, high-density liquid-like (HDL) and the more tetrahedrally ordered, low-density liquid-like (LDL) regions. These minima serve as integration boundaries for estimating the areas associated with each population.

In LiCl solutions, an additional peak is observed at low $\xi$ values ($\xi\lesssim 0.2$), corresponding to highly disordered local structures. This feature is excluded from consideration in the HDL and LDL population analysis. In such cases, the HDL region is defined between the first and second minima, while the LDL region lies beyond the second minimum. The fractions of HDL and LDL structures are obtained by integrating their respective regions and normalizing by the combined area.

Appendix F NEMD Simulation Results

The results of the NEMD simulations are given in Table LABEL:table:NEMD-data.

Appendix G Equilibrium NPT Simulation Results

The Equilibrium NPT simulation results are presented in Table LABEL:table:NPT-data.