Introduction

Worldwide water scarcity is a severe problem that has worsened in recent decades1. The factors exacerbating water scarcity include population growth, inefficient use of water resources, increased water usage in agriculture and industry2 and climate change3. Desalination technologies could mitigate water scarcity so have been receiving significant attention since the 1960s. In 2017, the global cumulative desalination capacity exceeded 100 million m3 day−11,4. To date, all desalination methods fall into two categories: materials-based methods such as reverse osmosis (RO), electrodialysis (ED)5 and ion-adsorbing materials6; and thermal-based methods such as multi-stage flash, interfacial evaporation7 and freeze desalination8. Thermal methods are appealing because they are driven by low-grade thermal energy (e.g. from waste heat or solar irradiation) and hence have the potential of decentralising desalination processes. However, current thermal desalination methods suffer from scale deposition, fouling and corrosion caused by phase change9, as well as high energy consumption at around 100 kWhth m−310. Preferred for their lower energy consumption between 3 to 7 kWhe m−311, RO has the highest installed capacity. ED recently started being developed at an industrial scale12,13 and can be more energy-efficient than RO when the salt concentration in the feedwater is lower than ca. 5000 ppm10, which is however much lower than the salinity of seawater (between 30,000 and 35,000 ppm). Despite being mature technologies, membrane fouling and degradation are still major issues related to all membrane methods14. Therefore, research interest in developing novel desalination concepts is growing4, e.g. capacitive deionisation15, ion-concentration polarisation16 and contactless steam generation17. Methods beyond the scope of desalination, e.g. atmospheric water harvesting18, are also receiving growing attention. Although most emerging technologies have not reached the maturity required for commercialisation because of their low throughput and high materials cost, their advantages such as decentralised desalination and use of low-grade thermal energy have justified further development.

Here, we re-visit a 150-year-old problem19, i.e. can thermodiffusion be used as an effective means for desalination? Thermodiffusion refers to species migration due to temperature gradients and is also termed the Soret effect or thermophoresis for large particles. Thermodiffusion is ubiquitous in nature and is thought to be related to many important phenomena including possibly the origin of life20,21. However, it has never been used as means for desalination, despite being a phenomenon first described in detail by Soret in 1879 through the observation of concentration inhomogeneities in a non-isothermal aqueous solution of NaCl19. There are only two theoretical papers briefly suggesting the possibility of thermodiffusion-based desalination22,23. Nonetheless, thermodiffusion has found use in few, yet essential, engineering applications, such as uranium enrichment in the Manhattan Project24, prediction of hydrocarbon distribution in oil reservoirs25,26 and analysis of biomolecular interactions27. In addition, many applications have been proposed based on thermodiffusion including nanoscale light tweezers28, carbon capture29, hydrogen separation30 and microfluidics-based separation of binary mixtures31,32,33,34. However, high-throughput applications based on thermodiffusion, such as thermodiffusive desalination (TDD), are yet to be achieved.

Thermodiffusion can be understood by considering the mass flux equation. In the presence of both concentration and temperature gradients, the mass flux in a binary solution is given by

$${{{{{{{\bf{J}}}}}}}}=-\rho D\nabla C-\rho C(1-C){D}_{{{{{{{{\rm{T}}}}}}}}}\nabla T,$$
(1)

where the first term on the right-hand side of the equation describes the isothermal (or Fickian) mass diffusion due to a concentration gradient C, and the second term describes thermodiffusion or the transport of species due to a temperature gradient T. In Eq. (1), D is the Fickian diffusion coefficient which is always positive because the molecules diffuse spontaneously from high to low concentration in a quasi-isothermal condition due to Brownian motion. In contrast, DT is the thermodiffusion coefficient which can be either positive (thermophobic) or negative (thermophilic), depending on thermal preference of the species35. Seawater is a multi-ion aqueous solution, but a quasi-binary approximation could be relatively accurate due to electroneutrality36. The Soret coefficient ST is defined as the ratio between the thermodiffusion and Fickian diffusion coefficients, ST ≡ DT/D. Based on Eq. (1), ΔC can be approximated as ΔC ≈ − C0(1 − C0)STΔT for quasi-linear temperature and concentration profiles with small changes. ST of NaCl in water is in the order of 10−3 K−137,38, meaning the ΔC that can be induced by thermodiffusion is small, limited by the ΔT achievable within the liquid phase of seawater. Based on measured values of ST, a temperature difference of 40 K yields ΔCwalls = 2400 ppm. between the top and bottom boundaries for an initial concentration of C0 = 30,000 ppm, i.e. 8% of C0. Therefore, thermodiffusion of NaCl in water can barely overcome Fickian diffusion, and this is likely why thermodiffusion has not been implemented as the separation principle in high-throughput applications. However, TDD is promising in that it is a moderate-temperature process that can be driven by temperatures much less than 100 °C extractable from the surrounding environment, e.g. from waste industrial heat or solar thermal energy. The heat being dissipated by cooling towers or sea currents could be used for desalination. Moreover, TDD is a completely materials-free separation process that relies on a simple concept. Note that thermodiffusion is stronger at higher concentrations with a peak efficiency when the mass fraction of the species is 0.5, as per thermodiffusion mass flux in Eq. (1). This suggests thermodiffusion could be a promising brine treatment method or precursor to existing materials-based desalination methods such as RO, ED, and emerging technologies such as ion-adsorbing materials6, which are all known to be more efficient with a lower-salinity feedwater.

In this work, we explore thermodiffusion as means of an effective water desalination method and reflect on its challenges and prospects. The method is termed thermodiffusive desalination, or TDD. First, we provide experimental results showing that thermodiffusion is able to achieve a tangible desalination level with a salinity reduction of 450 ppm and 700 ppm from an initial NaCl concentration of 30,000 ppm (seawater) and 60,000 ppm (brine), respectively. Importantly, the experimental setup is a simple parallel-plate channel with laminar flow, yet the concentration drop is relatively small after a single pass. As a proof of concept for TDD, re-circulation (multi-pass flow) through the same channel is implemented to increase the concentration drop to more than 2000 ppm. The multi-pass experiment was also performed using a substitute of seawater (i.e. a multi-ion aqueous solution) in which thermodiffusive separation of the cations was experimentally observed and then corroborated via molecular dynamics simulations. Finally, we propose and theoretically assess a cascaded structure, i.e. single flow pass through a multi-channel device, as a solution for scaling up TDD. We show a theoretical concentration drop exceeding 25,000 ppm with a recovery rate of 10%, which is amenable for high-throughput desalination. While still a proof of concept, this study provides valid pathways towards realising a single-phase thermal desalination process, demonstrating that water desalination is possible without the need for materials or phase change. Current materials-based methods utilise functional materials, such as ion-adsorbents6, selective membranes39 and permeable solar absorbers7, but these methods have intrinsic issues with fouling and regeneration. Avoiding the phase changes inherent in current thermal desalination methods reduces the relatively large latent heat requirements, giving TDD the potential to be less energy intensive without the cumbersome maintenance of functional materials.

Results

Thermodiffusive desalination concept

The proposed concept of TDD follows the configuration depicted in Fig. 1a. A multi-component salt solution passes through a rectangular channel with a vertical positive temperature gradient (heated from the top), as shown in Fig. 1b. This channel could be a simple, single channel—as in the thermodiffusive desalination unit (TDU) depicted in Fig. 1c—or a more complex multi-channel device (section “Burgers cascade”). Note that the species in this solution can be either thermophobic or thermophilic, depending on the local temperature and concentration. Salt ions in aqueous solutions are generally thermophobic above a temperature threshold called the inversion temperature T0 (ca. 12 °C for NaCl)37. In the designed TDU (Fig. 1c), a fully-developed laminar flow between thermostatically-controlled parallel plates is ensured, as transient flows above the critical Reynolds number induce mixing. A temperature difference is applied across the channel height h, with the top being heated and the bottom cooled to obtain a thermally stratified condition (otherwise, Rayleigh–Bénard instability could stir the flow40). In thermophobic transport, ionic species tend to accumulate towards the cold bottom side, resulting in a negative concentration gradient. Together with the negative density gradient brought by thermal stratification, species stratification ensures that natural convection-induced mixing does not occur in our thermodiffusive separation channel35. Heating can be achieved with low-grade heat including solar thermal energy harvested with solar absorbers41 or industrial waste heat. Cooling can be achieved through convection with the low-temperature saline water reservoir. The thermodiffusive separated solution is bifurcated into two streams at the end of the rectangular channel using a sharp spacer. In Fig. 1c, a more detailed view of the TDU with spacer (in the inset) is shown (see photos in Supplementary Fig. 1). Here, the volumetric flow rate Q of the saline feedwater is controlled with a peristaltic pump. The water is degassed before entering the TDU to avoid microbubble-induced mixing. Two hollow nickel-plated copper blocks are thermostatically controlled with a counter-flow water circulation system to establish a positive temperature gradient throughout the channel. Two bottles are placed at the same height so that equal flow into each bottle is ensured.

Fig. 1: Concept of thermodiffusive desalination and unit design.
figure 1

a Concept figure showing a laminar flow of saline water passing through a thermodiffusive separation channel. The temperature difference ΔT can be established with low-grade thermal energy. For temperatures above the inversion temperature, the thermophobic ions in the solution migrate towards the cold side as the saline flow progresses along the channel. This results in an upper stream having a lower salinity than the feedwater. The hydrodynamic and thermal entrance lengths, Lv and LT, are both less than 1% of the total channel length, Lcha. b The Reynolds number, Ra, indicates that Lv and LT are negligible compared to Lcha. Thus, the flow in the thermodiffusive separation channel can be approximated as a laminar, fully-developed planar Poiseuille flow with a positive quasi-linear temperature profile. x is the saline water flow direction and x = 0 corresponds to the channel inlet. The Péclet number, Pe, indicates an advection-dominant mass transport in the x direction and diffusion-dominant transport in the y direction. At the inlet, the solution is homogeneous, whereas at the outlet it becomes heterogeneous due to thermodiffusion. c Thermodiffusive desalination unit (TDU) design. The volumetric flow rate Q of the saline mixture in the TDU is controlled by a peristaltic pump between 1 to 16 mL min−1. Feedwater is degassed before entering the channel. The fluid path for the saline water is indicated by yellow lines. The channel is 500 mm long, 20 mm wide, and 1 mm high. At the exit of the channel, the saline water is separated into two streams by a spacer, as shown in the inset. After bifurcation, the two streams are collected in bottles at the same height, ensuring equal pressure and hence equal flow rate. Two hollow copper blocks with water circulation create the ΔT. Thermocouples were embedded in top and bottom channel walls to monitor the temperature (Supplementary Fig. 7).

Single-pass TDD

Since thermodiffusion in multi-component solutions is not adequately understood and there are limitations in the accurate measurement of diffusion in multi-component mixtures42,43, we started our investigation with a single-pass TDD experiment using the most common binary approximation of seawater: a NaCl/H2O solution. We recently reported numerical simulation results22 of the thermodiffusive separation within a parallel-plate channel with a plane Poiseuille flow and a linear temperature gradient. Details of the numerical simulations relevant to this work are available in the Supplementary Method 1, along with Supplementary Figs. 25. The design rationale for the TDU is detailed in Supplementary Method 2, along with Supplementary Fig. 6. The flow within the TDU is laminar (Re < 350) because of the reduced flow speed required to achieve a nearly fully developed concentration profile at the channel outlet, i.e. a rather long resilience time (t > 1 min) within our 0.5 m channel is needed to obtain a nearly complete thermodiffusive separation. With our TDU design (h = 1mm) and _target flow rates (1 – 16 mL min–1), the hydrodynamic, thermal, and concentration entrance lengths are around 0.1%, 1%, and 100% of the entire channel length, respectively. We also found that the mass transport is advection-dominant along the channel (Péclet number Pex ≈ 103) while diffusion and thermodiffusion dominate across the channel height (Pey = 0).

The experimental and numerical results are shown in Fig. 2a, where the concentration difference ΔC (between top and bottom solutions in the collection bottles of Fig. 1c) is shown as a function of the measured temperature difference ΔTmeas along the channel (measured with thermocouples). A typical temperature profile along the channel is shown in Supplementary Fig. 7b. The mean temperature Tmean and temperature difference between walls ΔT were chosen such that the temperature of the bottom wall Tcold in the TDU (either single channel or a Burgers cascade) is above the inversion temperature at which thermophobic–thermophilic transition occurs. The inversion temperature has been reported to be between 10 °C38 and 12 °C37 for an aqueous NaCl/H2O at seawater concentration. Measurements and modelled predictions of ΔC were conducted via phase-shifting interferometry (PSI)44 and computational fluid dynamics (CFD)22, respectively. An accurate measurement of the ΔC is critical in assessing the performance of TDD, and we found that our in-house PSI was more accurate than commercial salinity meters based on electrical resistance. Supplementary Method 3 and Supplementary Fig. 8 provide details on the PSI measurement and associated uncertainty. The vertical error bars in Figs. 2a, b and 3a, b depict the errors δΔC in concentration difference ΔC measurements, and the errors δCdrop in concentration drop Cdrop measurements, respectively. Measurement errors mainly arise when extracting the phase difference from unwrapped phase-shifted data (insets of Fig. 2a) between Chigh and Clow solutions (i.e. those extracted from the collection bottles in Fig. 1c). A typical relative PSI-based error \(\frac{\delta \Delta C}{\Delta C}\) is around 10% for the relatively small salinity difference ΔC < 1000 ppm. In contrast, a commercial electrical resistance-based salinity meter may have a relative error exceeding 100% since the measured ΔC is very small. The CFD model takes the measured temperature along the top and bottom walls of the channel (Supplementary Fig. 7) as boundary conditions and considers temperature-dependent Soret and Fickian diffusion coefficients, while assuming a fully developed flow between the plates.

Fig. 2: Performance assessment of single-pass thermodiffusive desalination.
figure 2

a The concentration difference ΔC between the top and bottom collection bottles as a function of the average measured temperature difference along the channel ΔTmeas for different mean temperatures Tmean. The experimental results (solid markers) are plotted with horizontal and vertical error bars representing the standard deviation of ΔT measurements (with six thermocouple pairs) and the uncertainty in ΔC measurements (with phase-shifting interferometry, PSI). The insets are PSI images following the experimental procedure reported in35,44. CFD results (dashed lines) with experimental temperature values as boundary conditions are included. Monotonic dependencies are revealed. b The effect of flow rate Q on ΔC is reported for Tmean ≈ 41 °C and ΔTmeas ≈ 34 K both experimentally and numerically. The errors in the experimental ΔC are standard deviations from PSI measurements. CFD Case 1 uses the thermocouple temperatures measured in each experiment as boundary conditions, whereas CFD Case 2 uses the common wall temperature profile of the experiment with lowest volumetric flow rate of 1 mL min−1. ΔC remains roughly the same when Q < 5 mL min−1. The effect of slight variation in the temperature field on the separation is prominent for small Q. c CFD Case 1 concentration contours for different ΔTmeas and Q. Larger ΔC are observed for greater ΔTmeas and smaller Q. d CFD computation of the concentration difference between collection bottles. Temperature-dependent ST yields a quasi-linear concentration profile within the cell (d1), together with the parabolic velocity profile (d2), yields the ion mass flux as a function of the height (d3). Dividing the integrated ion mass flux in the top and bottom halves with the corresponding flow rate yields the concentration at the collection bottles.

Fig. 3: Multi-pass thermodiffusive desalination for aqueous NaCl and seawater.
figure 3

Both seawater level and brine level (double the seawater concentration) are considered. a The initial concentration of the NaCl in aqueous solution was 60,000 ppm. NaCl concentration in both the top and bottom streams are plotted for each pass. Both experimental results and numerical results are demonstrated. The errors in the experimental ΔC are standard deviations from PSI measurements. In each pass, the top stream with a lower salinity was accumulated and put into the next pass. The NaCl concentration in the top stream kept dropping, demostrating that thermodiffusion is scalable beyond a single pass. b The concentration drop is plotted as a function of number of passes recirculated for two different initial concentrations (30,000 and 60,000 ppm). On the left axis, the concentration drop Cdrop achieved after four passes were 1200 ppm and 2000 ppm for an initial concentration of 30,000 and 60,000 ppm, respectively. On the right axis, the recovery rate is shown. There is a trade-off between Cdrop and recovery rate. c Same experiment was performed using multi-ion seawater substitute. The concentration of different cations: Na+ (c1), Mg2+ (c2), Ca2+ (c3), and K+ (c4), was measured by two different ICP methods for the same set of samples collected in each pass. Linear fits were performed for the top stream concentrations to show the concentration reduction trend. The errors in ICP-AES measurements are the standard deviations between five replicates. The errors of the ICP-OES results were not available. Na+ is most abundant in quantity and provides most reliable results. For ICP-AES results, the Ctop = (−190n + 14,200) ppm, i.e. Cdrop = 190 ppm and C0 = 14,200 ppm. For ICP-MS results, Ctop = (−111n + 13,329) ppm. Linear fits are also applied to Mg2+ and Ca2+. For K+, however, no trendline is fitted due to lack of discernible trend.

For the same mean temperature Tmean, we could confirm that ΔTmeas is roughly proportional to ΔC, which is an expected result according to Eq. (1). Insets of Fig. 2a show clear fringe patterns between the two outlet samples, which is evidence of thermodiffusive separation of NaCl in the channel. Under a larger ΔTmeas, more fringes in the phase-shifted data are indicative of a larger ΔC between the top and bottom solutions. Furthermore, the good agreement between PSI and CFD results confirm that no mixing occurred in our channel (as the latter assumes a fully-developed laminar flow). Note that, for the same ΔTmeas, ΔC increases with mean temperature Tmean. This is expected as the Soret coefficient of aqueous NaCl is known to increase monotonically with temperature37.

The effect of volumetric flow rate Q was also analysed. Figure 2b shows that a near-maximum ΔC was achieved for Q < 5 mL min−1. This is because the saline water from the inlet remains in the 1 mm high channel for more than 2 min, which is enough time for full separation to occur. Larger flow rates result in a drop of ΔC as the residence time of the solution in the channel is not long enough for thermodiffusion to fully separate NaCl before reaching the spacer. Furthermore, different water bath flow configurations yield different wall temperature profiles (Supplementary Fig. 7). A counter flow configuration with the hot water bath flowing in the opposite direction to the saline water was chosen (Supplementary Method 4). Bubbles from the water circulation may accumulates at the heat exchange surface (Supplementary Fig. 1b), which reduce the local heat transfer coefficient and hence lower ΔT. Without visible bubbles, the ΔC obtained for the same set water circulation temperatures is within the PSI measurement accuracy, meaning an excellent repeatability was obtained. However, despite aiming at having the same wall temperature profiles when investigating the effect of the volumetric flow rate Q, the wall temperature varied due to the different thermal resistance caused by the bubbles. Difference in wall temperature slightly change predictions of ΔC for CFD Case 1 and Case 2, as shown in Fig. 2b when the volumetric flow rate is less than 5 mL min−1. CFD Case 1 was performed using measured wall temperature profiles (for the corresponding Q) as boundary condition, while CFD Case 2 was performed using the wall temperature profile obtained in the Q = 1 mL min−1 experiment. The impact of slight variations in temperature profiles on the ΔC is more obvious for small Q. In Fig. 2c (CFD Case 1), concentration profiles on the vertical plane along the channel are shown for three sets of Q and ΔTmeas. As expected, a smaller separation is obtained for a smaller ΔTmeas (Fig. 2c1 vs. Fig. 2c2, c3). Furthermore, when a very small Q = 1 mL min−1 is implemented (Fig. 2c2), the concentration profile becomes fully developed halfway through the channel (slight variations are still seen in Fig. 2c2, but these are due to non-uniform wall temperatures as shown in Supplementary Fig. 7). In contrast, a larger Q increases the concentration developing length (Fig. 2c3).

The calculation procedure for obtaining the concentration at the top and bottom collection bottles from CFD results is depicted in Fig. 2d. We show that despite the seemingly large ΔCwalls between the top and bottom boundaries, due to the parabolic velocity profile, the ΔC between the top and bottom collection bottles is significantly smaller than in a convectionless environment. The width of the TDU channel is w and the top stream concentration is calculated as \(\frac{\rho w}{Q/2}\int\nolimits_{h/2}^{h}C(y)u(y){{{{{{{\rm{d}}}}}}}}y\), whereas the bottom stream concentration as \(\frac{\rho w}{Q/2}\int\nolimits_{0}^{h/2}C(y)u(y){{{{{{{\rm{d}}}}}}}}y\). The ΔC between top and bottom streams is around 37% of that between the boundaries. Hence, for ΔC = 900 ppm as in Fig. 2a, the ΔCwalls between top and bottom boundaries is around 2400 ppm.

Multi-pass TDD

The maximum ΔC after a single pass (Fig. 2a) is ca. 900 ppm for ΔTmeas = 37 K. However, 900 ppm is only 3% (relative value) of the initial NaCl concentration in seawater, ca. 30,000 ppm. Moreover, the salinity reduction or concentration drop Cdrop of the top stream concentration from the inlet value is only half of ΔC between top and bottom streams (i.e. 450 ppm). This small drop is not enough to obtain low-salinity water useful for agriculture, which may require 95% relative salinity reduction and accounts for 69% of water use worldwide45. To amplify thermodiffusive separation, a simple way is to accumulate the top low-concentration stream in a container and then pass it again through the same channel. This re-circulation experiment was done for two initial NaCl feedwater concentrations, seawater level (30,000 ppm) and brine level (60,000 ppm). The concentration of NaCl in the top and the bottom streams after each pass is shown in Fig. 3a for brine feedwater. Note that the bottom high-concentration stream is discarded after each run, while the top low-concentration stream is re-circulated. CFD Case 1 results have a very good agreement with the experiments, except for the last pass as the recovered volume was low and contamination could have affected the experimental result.

Figure 3b shows the concentration drop Cdrop increasing with the number of passes, demonstrating that TDD can be scaled up beyond the modest results obtained from a single pass. Furthermore, it is noticed that Cdrop is larger in brine than in seawater, which can be attributed to the larger thermodiffusive mass flux at increased concentration following Eq. (1). This feature suggests that TDD may be a potential pre-treatment method for high-salinity feedwater. However, the significant drawback of this multi-pass TDD is that the volume of the desalinated stream drops exponentially, halving after each pass. In the fifth pass, the recovery rate was ca. 4%, which is excessively small for a yield with relatively large salinity, ca. 95% of its initial value. Therefore, another method for amplifying thermodiffusive separation should be devised for it to be amenable in desalination applications.

Another factor that may affect TDD in real applications, is that seawater is a multi-component ionic solution and some ions may hinder thermodiffusive separation if their thermodiffusive transport opposes or is weaker than that of the thermophobic ions of Na+ and Cl (the major species in seawater). Despite being discovered more than 150 years ago, the mechanisms behind thermodiffusion are still unclear and quantitative measurements have been limited to ternary mixtures46. Here, we are more interested in the ‘collective’ thermodiffusive transport of all the seawater components assuming that electrostatic interactions have them diffuse as a single effective species. To assess TDD in seawater, natural sea salt was dissolved in deionised water at a concentration of 30 g L−1 to prepare an artificial substitute of seawater. Figure 3c shows the multi-pass experimental result with seawater substitute. Inductively coupled plasma atomic emission spectroscopy (ICP-AES) and inductively coupled plasma mass spectrometry (ICP-MS) were used to measure the concentration of different cations in the same set of samples collected. Although the exact composition of the seawater substitute was not the same as in natural seawater reported in the literature10, it demonstrates that the TDD principle applies to multi-ion aqueous solutions.

A linear fit was applied to the top solution concentrations to determine the concentration drop Cdrop in each pass. Based on the numerical analysis in Fig. 2d, Cdrop can be translated back to the concentration profile across the channel height to calculate the Soret coefficient ST. The process to derive ST from Cdrop is detailed in the Supplementary Method 5. ST was \({2}_{-1}^{+1}\times \,1{0}^{-3}\,{{{{{{{{\rm{K}}}}}}}}}^{-1}\) based on both the ICP-AES and ICP-MS measurement results. The superscripts and the subscripts originates from the process of performing a linear fit to Cdrop(n) (Fig. 3c). For the linear fit, the slope is the concentration drop per pass Cdrop and the intersection at n = 0 is the initial concentration C0. When assuming a planar Poiseuille flow (Fig. 2d2), both ΔCwalls and the concentration profile C(y) (Fig. 2d1) can be calculated from Cdrop (Fig. 2d3). With the known C(y), ST can be calculated as per Supplementary Method 5. The superscripts and subscripts correspond to the linear fit with the largest and smallest slopes with a confidence interval of 95%, respectively. These values are comparable to or greater than that of the binary NaCl/H2O. Although an accurate measurement of ST from ICP analyses is challenging, we can conclude that TDD is not hindered by the multi-component nature of seawater.

TDD at the molecular level

Molecular dynamics (MD) simulations were performed on NaCl and multi-ion seawater substitute solutions to assess the application of the TDD principle to multi-component electrolytes. The TDU experimentally achieved a salinity reduction in both NaCl (Fig. 3a, b) and seawater substitute solutions (Fig. 3c). However, there was a large uncertainty in the reported salinity measurements using ICP methods. Therefore, simulations at the molecular level are sought to confirm the experimental observation that thermodiffusive separation occurs in a multi-component ionic solution of seawater. Furthermore, MD simulations may also provide insight into ion–ion and ion–water interactions, as well as providing a single theoretical framework to compare thermodiffusion in artificial seawater substitute and the simpler binary NaCl/H2O approximation.

The MD simulation results shown in Fig. 4 predict seawater ions exhibit a thermophobic behaviour (see Supplementary Method 6 and Supplementary Figs. 911 for more details on methods and MD results). The simulation was run with an average temperature of 40 °C and a temperature difference of 40 K to match the optimal operating conditions of the TDU and maximise the expected concentration drop (see Fig. 2a). The total concentration of the multi-ion seawater model and NaCl solutions was more than doubled (i.e. brine) to improve the counting statistics of the MD results. The exact compositions of the seawater brine is available (Supplementary Table 1). The tested model solutions are therefore referred to as seawater brine and NaCl brine. The average concentration of seawater ions shown in Fig. 4a increases as the temperature in the simulation decreases, demonstrating that the thermodiffusion of seawater ions is thermophobic in the temperature range 20–60 °C.

Fig. 4: Molecular dynamics modelling of thermodiffusive desalination.
figure 4

a Molecular dynamics (MD) model of seawater brine under a non-uniform temperature profile spanning ca. 20–60 °C. The ion concentration (for all ions) decreases as the temperature increases, demonstrating thermophobic behaviour of seawater ions in silico. The simulation volume (top) contains water molecules (light grey), and seawater ions (multiple colours). A cold and hot thermostat of width 1 nm are placed at x = 5 nm and x = 15 nm (highlighted regions). The steady state temperature (solid line) and time-averaged ion concentration (markers) between \(\left[{t}_{0}=20,\,{t}_{{{{{{{{\rm{final}}}}}}}}}=280\right]\) ns of simulation time are plotted. The error bars are the standard deviations for the time-averaged temperature and concentration within the steady state simulation time. b Average ion concentration in simulated seawater brine and NaCl brine across temperature. The concentration drop for seawater brine (7000 ± 400 ppm) was larger than concentration drop for NaCl brine (3800 ±  600 ppm). The ion concentration for seawater brine (in orange) and NaCl brine (in blue) is averaged over four and three independent replicates, respectively, and shown with standard error in the mean of the local temperature (markers with error bars). The spread in the average ion concentration is shown at a 95% confidence interval (shaded regions) due to the small number of replicates. A linear line of best fit is reported with an R2 value for seawater brine (solid line) as C(T) = (−176 ± 11)T + (78,600 ± 500) ppm and NaCl brine (dashed line) as C(T) = (−95 ± 14)T + (73,600 ± 600) ppm. c Comparison of Soret coefficients obtained from MD simulations and experiments (both from this work and literature). (i) The Soret coefficient calculated for seawater brine is 1.8 times larger than that of NaCl brine. The errors arise from the ST calculation process. When applying different possible fits to the results (with a confidence interval of 95%), they result in different ΔC. Details of the error derivation are available in Supplementary Method 5. The double asterisk annotation (**) indicates that the increase in the Soret coefficient is statistically significant with p = 0.003 (i.e. p < 0.01), reported by a one-sided t-test. (ii) The Soret coefficient for NaCl brine obtained from MD simulations is comparable to that from experiments. (iii) The experimental Soret coefficient for NaCl in this work is lower than that of ref. 37.

MD simulations predict that the modelled multi-component seawater had a greater salinity reduction than that in the binary NaCl/H2O system. Both brine of seawater and NaCl behaved thermophobically, with ion concentration decreasing with increasing temperature, as shown in Fig. 4b. Seawater brine had a ΔCwalls of 7000 ±  400 ppm between the hot and cold boundaries, larger than the 3800 ±  600 ppm ΔCwalls predicted for NaCl brine. A comparison of the Soret coefficients in Fig. 4c(i) shows that the predicted Soret coefficient of \(2.{4}_{-0.1}^{+0.0}\times \,1{0}^{-3}{{{{{{{{\rm{K}}}}}}}}}^{-1}\) for seawater brine was \(1.{8}_{-0.2}^{+0.3}\) times larger than the predicted Soret coefficient of \(1.{4}_{-0.2}^{+0.2}\times \,1{0}^{-3}{{{{{{{{\rm{K}}}}}}}}}^{-1}\) for NaCl brine. To relate the scatter plot as in Fig. 4b with ST, a linear fit was performed to obtain a continuous concentration profile C(y). ST is then calculated by finding the value of ST that renders an analytical solution of C(y) close to the linear fit. The linear fit with the largest and smallest slopes corresponds to the superscript and subscript in the reported Soret coefficient, respectively. Details of the analytical process can be found in Supplementary Method 6. Analysis of the individual ion species concentration profiles suggest that the larger separation predicted for seawater brine is caused by the stronger thermodiffusive separation of Mg2+ and the consequent electrostatic interaction between Mg2+ and the other ions (Supplementary Fig. 11). Overall, the MD simulations predicted that thermodiffusion can produce a greater salinity reduction in our multi-component seawater brine model than in NaCl brine.

Confidence in the enhanced thermodiffusive response of multi-component seawater brine is gained from the agreement between MD simulation predictions and experimental measurements of NaCl brine. The 3800 ± 600 ppm ΔCwalls of NaCl brine predicted by simulation (reported in Fig. 4b) is larger than the 3100 ± 160 ppm ΔCwalls of NaCl brine measured experimentally in the TDU (reported in Supplementary Table 2) for comparable concentrations. However, note that the MD simulation has a temperature difference of 40 K, and the comparable NaCl experiment achieved a reduced temperature difference of ca. 31 K (Supplementary Table 2). At this reduced temperature difference, the MD result predicts a ΔCwalls of 2900 ± 400 ppm (calculated from the best linear fit for NaCl brine in Fig. 4b), better agreeing with the TDU-measured ΔCwalls. The Soret coefficient of NaCl brine predicted by simulation is \(1.{4}_{-0.1}^{+0.2}\times \,1{0}^{-3}{{{{{{{{\rm{K}}}}}}}}}^{-1}\), a slight underestimate compared to the experimentally determined Soret coefficient of \(1.6{5}_{-0.14}^{+0.08}\times \,1{0}^{-3}{{{{{{{{\rm{K}}}}}}}}}^{-1}\), as reported in Fig. 4c(ii). However, the estimates are within their error ranges. Experimental values for the Soret coefficient found in this work are lower than that of the literature value for the Soret coefficient published in37, reported in Fig. 4c(iii). Overall, the simulation predictions for NaCl brine agree with experimental measures of salinity reduction and Soret coefficient made in this work, and are less compared to early literature values37. As there are no comparable values for the Soret coefficient of the multi-ion seawater substitute published in the literature, validation of the NaCl brine Soret coefficient gives confidence to the enhanced thermodiffusive response predicted for seawater brine.

Burgers cascade

The simple multi-pass re-circulation scheme (Fig. 2) does not meet requirements of practical desalination applications as it has a very low recovery rate. Such configuration is not efficient as it discards solution that is already of a lower concentration than the initial feedwater. Here, we further develop TDD with a scalable thermodiffusive separation device known as the Burgers cascade47, previously shown to enhance thermodiffusive separation in binary gas mixtures30 but never in liquid systems. The working principle is shown in Fig. 5a. The Burgers cascade is a device with multiple top(Clow)/bottom(Chigh) vertical bifurcation and right–left horizontal recombination cells, each essentially acting as a miniature TDU. The possible dimensions and design of the cell is available in the Supplementary Method 7. M and N refer to the total number of rows and columns in the Burgers cascade, respectively, while m and n are the corresponding indices. This complex interconnection of many small channels may be blocked by two-phase formation. However, we demonstrated that bubbles can be avoided by slightly tilting the Burgers cascade and filling it up from the bottom. A miniature Burgers cascade with M = 5 and N = 10 was manufactured and the flow structure was visualised with dye as shown in Fig. 5a2. More details of this flow visualisation miniature is available (Supplementary Fig. 12).

Fig. 5: A scalable cascaded TDD device.
figure 5

a The structure of the Burgers cascade. a1 The Burgers cascade contains many small cells; each cell is a small parallel-plate channel. In a single cell, the working fluid is divided into two streams, top and bottom, denoted by red and blue arrows respectively. For the entire Burgers cascade, the top is hot and the bottom is cold. The hot stream (red arrow) flows towards the top-left cell outlet and has a lower concentration for thermophobic species, while the cold stream (blue arrow) flows to the bottom-right outlet and has a higher concentration. Water flows in the m direction, where m is a row and n is a column. a2 Photo of a transparent miniature with M = 5 and N = 10 showing the lateral spread of dye after injection close to the device inlet. Details of this experiment are provided in Supplementary Fig. 12. b Contour plot of the concentration in each cell. C0 = 30,000 ppm, ΔT = 60 K, ST is assumed to be 1.8 times of that of NaCl in water. M = 185 and N = 20 (b1) and M = 490 and N = 35 (b2). When the cut-off distance is 10% of the total distance, i.e. the recovery rate Rw is 10%, the yield concentration Clow = 5000 ppm (b1) and 1000 ppm (b2), respectively. c The concentration of the cold stream is plotted along the flow direction (when N = 20) for different recovery rates. Two different feedwater concentrations are considered, C0 = 30,000 and 60,000 ppm, with the same assumptions as in (b) were applied. Drinking water here is defined as water with salt concentration less than 1000 ppm. Lower recovery rate means larger ΔC. d The yield concentration when expanding the size of the Burgers cascade fixing Rw to 10%. By expanding the size of the Burgers cascade, seawater can be desalinated to drinking water standard. For different levels of _target salinity, the corresponding size M is indicated by coloured labels.

We assume that seawater behaves as a binary mixture of ions and water. Based on the MD work, seawater ions are predicted to have an effective ST 1.8 times larger than that of NaCl/H2O solution. Other assumptions are: mixing of the two streams happens instantaneously at the inlet of each cell; the concentration profile fully develops in each cell; the two streams at the outlet of the cell are perfectly bifurcated at the horizontal mid-plane. After applying our numerical model for a single-pass binary solution to each cell, we can obtain the concentration in each cell of the Burgers cascade. Figure 5b shows the contour plot of concentration in each cell in the Burgers cascade when ΔT = 60 K is applied. The recovery rate can be adjusted by selecting the cut-off location at the Burgers cascade outlet. With 20 cells in each row (N = 20), Clow of the yield stream when varying M is plotted in Fig. 5c. With 20 cells in each row and 185 rows, which is comparable to the Burgers cascade with 22 cells in each row and 450 rows in a previously demonstrated experiment30, the salinity reduction is 50% even when the recovery rate is as high as 50%, for both seawater and brine concentrations.

The Burgers cascade can also be scaled to produce different yield concentrations, as shown in Fig. 5d, with the required number of cells indicated in the legend. Following an analytical method (Supplementary Methods 7 and 8), the pressure drop in the Burgers cascade can be calculated. Even for the case when M = 490 and N = 35 (Supplementary Fig. 13a), the pressure drop is much smaller than 1 bar, which is what a modest pump can supply and less than 5% of RO pumping requirements48. It is important to note that for most designs listed, the pressure drop is less than 10 kPa and the electrical energy consumption to produce 1 m3 of fresh water, when considering only the pump work to overcome the pressure drop, is less than 3 Whe (Supplementary Method 8). The minimum thermodynamic limit for the equivalent case49, i.e. with a feedwater concentration of 30,000 ppm, yield concentration of 5000 ppm and recovery rate of 10%, is 0.4 kWh m−3 (Supplementary Method 8). Therefore, the electrical energy consumption of TDD alone is significantly lower than the theoretical minimum energy of separation, only 1% of its value for the aforementioned condition. Moreover, the mass flux equation Eq. (1) dictates that thermodiffusive mass flux monotonically increases with concentration until C = 0.5. Combined with the low electrical power consumption of the Burgers cascade, TDD has the potential of being used in a hybrid configuration together with other desalination technologies that are more efficient when the feedwater salinity is low, e.g. RO and ED. When used in conjunction with RO, for example, considering only the electrical energy consumption, the overall energy consumption is reduced by more than 80% from 4.5 kWhe m−3 to 0.7 kWhe m−3 (Supplementary Method 8). TDD may also serve as a pre-treatment method for ED as it is more energy efficient than RO when the feedwater has a salinity lower than 5000 ppm10,50. Last but not least, TDD by itself has a comparable salt removal rate as compared to other emerging desalination methods (Supplementary Table 3), justifying its further development beyond the proof of concept presented in this paper.

TDD has interesting trade-offs, such as that between the flow rate Q and the thermal energy consumption (Supplementary Fig. 13b). The concentration profile development time is inversely proportional to the cell height squared, i.e. smaller cell height is preferred to obtain a larger yield (Supplementary Fig. 5c). On the other hand, a larger cell height is preferred to minimise heat flux through the structure because when ΔT is fixed, the heat flux across the structure is inversely proportional to the channel height (Supplementary Method 7).

Discussion

Proof of concept

This work presents a theoretical and experimental framework for implementing thermodiffusion in single-phase thermal desalination applications. We provide modelling results and lab-based measurements that show thermodiffusion can reduce the salinity of seawater. The salinity reduction in a single-pass TDD channel is 450 ppm for a temperature difference of ca. 37 K and mean temperature of 40 °C. Although the salinity drop for a single-pass TDD is small, note the relatively low temperature for operating TDD: this temperature range is ubiquitous as waste heat in industrial settings or hot deserts where water is needed. We experimentally demonstrated that TDD can be scaled up by re-circulating the low-concentration stream through the channel and achieved a salinity reduction (after four passes) of 1200 and 2000 ppm for an initial feedwater concentration of 30,000 ppm (seawater) and 60,000 ppm (brine), respectively. The experimental results are in agreement with our numerical modelling in the continuum regime.

In addition, the thermodiffusive separation in a multi-component ionic aqueous solution was quantified. Due to large errors in ICP spectroscopy, accurate concentration measurements for each ion were not attainable for the multi-component seawater substitute. Nonetheless, our measurements indicate that an effective Soret coefficient is likely to be larger in seawater than in binary NaCl/H2O solutions. Such result agrees with our MD simulation results that show the effective Soret coefficient is larger in multi-ion brine than NaCl brine. The salinity reduction of natural seawater via TDD is therefore expected to be larger than reported values for NaCl/H2O solutions.

We recognise that TDD in the simple parallel-plate channel structure does not seem feasible for practical applications requiring a salinity lower than 5000 ppm and a large recovery rate. As a potential solution to scaling up TDD, we proposed and analysed a multi-channel bifurcation and recombination device called Burgers cascade30. Despite the complex internal structure, the pressure drop is not a major concern as the flow rate is very small. We report a theoretical concentration drop greater than 25,000 ppm via thermodiffusive separation using the Burgers cascade with a recovery rate of 10%.

Advantages

TDD has two major benefits. First, TDD is a thermal desalination process entirely operated in the liquid phase without relying on any functional materials such as membranes or ion-adsorbents that need regular maintenance to avoid extensive materials degradation. It is worth noting that all the TDD experiments reported here were run intermittently for more than one year using the same TDU device (depicted in Supplementary Fig. 1). There was no visible corrosion because the copper blocks were plated in nickel to avoid such surface degradation. From a process perspective, TDD is most similar to ED, where Fickian diffusion mixing is overcome by applying a competing phenomenon, which in this case is thermodiffusion. Second, TDD is driven by relatively low-temperature heat that can be sourced from readily available sources including industrial waste, solar irradiation or simply the surrounding environment. In Fig. 6, a comparison of the specific energy consumption for different desalination technologies is presented. The references corresponding to each numbered item are listed in Supplementary Table 4. We see that thermal-driven methods usually have a much higher energy consumption. Nonetheless, a higher energy consumption can be justified if it is from a source that is readily available in the environment or if it is waste heat (e.g. from industrial or agricultural processes). The use of low-grade thermal energy as the driver for desalination has motivated the development of various desalination technologies including membrane distillation, novel solar-driven distillation and TDD (this work). From an energy perspective, TDD is most similar to membrane distillation where waste heat and large temperature gradients are often implemented. The electrical energy required to produce 1 m3 of 5000 ppm yield water (concentration drop of 25,000 ppm) is less than 3 Whe (Supplementary Method 8), which is significantly lower than the minimum thermodynamic limit for desalination under the same recovery and concentration conditions. This indicates that TDD has an outstanding potential in areas where surrounding thermal energy is abundant while access to electricity is limited.

Fig. 6: Comparison of different desalination technologies.
figure 6

The energy consumed per unit volume of yield, or specific energy consumption (SEC), is plotted against the concentration drop for different desalination technologies: reverse osmosis (RO), multi-stage flash (MSF) and multi-effect distillation (MED), electrodialysis (ED), adsorption, capacitive deionisation (CDI), novel solar-driven desalination (NSD), e.g. interfacial evaporation and contactless steam generation, membrane distillation (MD), and the thermodiffusive desalination (TDD) technology proposed in this work. Corresponding sources are indexed and listed in Supplementary Table 4. For some technologies, SEC is dependent on the feedwater salinity and the yield salinity. The vertical black dashed line indicates a concentration drop from seawater to freshwater. The thermal-driven desalination methods are captured by the mauve circle. SEC of different technologies vary vastly and the thermal-driven desalination technologies generally have a much higher energy consumption. For TDD, the presented SEC data is based on multi-pass TDU experiment. A heat flux of 232 W (Supplementary Method 8) was used and the yield was 0.3 L h−1. Under this condition, SEC is calculated as \({{{{{{{\rm{SEC}}}}}}}}=0.232\,{{{{{{{\rm{kWh}}}}}}}}/0.3\,{{{{{{{\rm{L}}}}}}}}=773\,{{{{{{{{\rm{kWh}}}}}}}}}\,{{{{{{{{\rm{m}}}}}}}}}^{-3}\). However, when only considering the electrical energy usage, \({{{{{{{{\rm{SEC}}}}}}}}}_{{{{{{{{\rm{e}}}}}}}}}\), its value drops to nearly 0. SEC data is not available for the adsorption-based desalination.

Challenges and future directions

There are two major challenges that limit the performance of TDD. First, there is a trade-off between yield and heat flux: a smaller height is preferred for a larger yield (Supplementary Fig. 6b) while a larger height is preferred for a smaller heat flux. Second, thermodiffusive separation has a fundamental limit indicated by the mass flux equation Eq. (1). The product C(1 − C) indicates that thermodiffusive mass flux diminishes with decreasing concentration, meaning that TDD is less effective when the feedwater concentration is lower. Considering these challenges, there are several pathways towards the practical implementation of TDD.

First, achieving a better understanding of the fundamental mechanism behind thermodiffusion is paramount to propose strategies that enhance thermodiffusive separation. For example, increasing the isothermal diffusion coefficient shortens the time for the concentration profile to develop in each cell, as the concentration development time is only related to the Fickian diffusion coefficient when the Soret coefficient is small (Supplementary Fig. 5a). Also, increasing the Soret coefficient would enable a greater salinity drop and thus a smaller number of cells is required in the Burgers cascade. Second, as described in Supplementary Method 8, TDD could be used as a pre-treatment method in a hybrid desalination process, concurrently implemented with more mature technologies such as RO and ED, which are more efficient at lower feedwater concentrations. Third, beyond the application of desalination, thermodiffusion can also contribute to the treatment of brine or industrial waste water that contains high-concentrations of pollutants. Thermodiffusive separation is in fact more efficient at salinity levels greater than seawater, as per Eq. (1). A further technology development for TDD would enable an accessible and environmentally friendly desalination method with comparable performance to other desalination concepts recently developed (Fig. 6 and Supplementary Table 3). While still a proof of concept, this paper presents theoretical and experimental underpinnings that may inspire further use of thermodiffusion as a separation principle in water treatment and other environmental engineering applications.

Methods

Simulation of thermodiffusive separation in a channel

Simulations were performed following the same procedure as in our previous work22. The channel is a simple parallel-plate channel. The continuum thermodiffusion modelling was based on the conservation of chemical species (NaCl) in a fluid flow with constant velocity field u following the governing equation: \(\frac{\partial }{\partial t}(\rho C)+\nabla \rho {{{{{{{\bf{u}}}}}}}}C=-\nabla {{{{{{{\bf{J}}}}}}}}\). The essential assumptions are: one-dimensional flow with a known velocity profile (a planar Poiseuille flow), linear temperature profiles across the channel height, and temperature-dependent ST and D37. The boundary conditions for temperature are matched to measured temperature values in experiments where twelve thermocouples (six on each side) are placed at the top and bottom boundaries of the 500 mm long channel. The temperature at each grid point is linearly interpolated from the thermocouple readings. The transient change in temperature readings is small and is ignored in the simulation. A fully-implicit finite volume method was used with the discretisation equations available in Supplementary Method 1.

Sample preparation

NaCl/H2O samples are prepared on a gravimetric basis using NaCl (Ajax Finechem, Thermo Fisher Scientific, minimum assay 99.7%) and deionised water (Pacific TII 12, Thermo Scientific). The balance has an accuracy of 1 mg (ATX224, Shimadzu). Artificial seawater was prepared using the same apparatus and procedure with rock sea salt sourced from Whyalla, South Australia (Rock sea salt, Olsson’s Salt), which was harvested through solar evaporation with no additives.

Thermodiffusive desalination experiment

In a single-pass experiment, a peristaltic pump (Kamoer FX-STP2) was used to control the volumetirc flow rate Q with an accuracy of 0.1 mL min−1 and a silicone membrane gas exchanger (PDMSXA-1000: PermSelect®1000) connected to −0.8 bar vacuum was used to degas the solution before it enters the channel. The copper blocks that accommodate water bath flow are nickel plated and corrosion was not observed before and after all the experimental runs. To monitor the temperature, six pairs of thermocouples were distributed 100 mm apart on both upper and lower plates of the channel. These K-type thermocouples (410-305, TC Direct) are read out by thermocouple readers (U6-Pro, LabJack). In multiple-pass experiment, between each pass, deionised water was used to flush the channel for around 15 min then compressed air was used to dry the channel for around 15 min to reduce the chances of contamination between solutions of different passes.

Optical measurement of concentration drop

Phase-shifting interferometry (PSI) is a digital interferometry technique that can resolve the phase difference ψ between test and reference beams. Our PSI technique uses a polarising Mach–Zehnder interferometer. The PSI layout was first developed in44,51. The laser (HNL050LB, Thorlabs, wavelength λ = 630 nm) is a polarised HeNe laser and a polarising beam splitter divides it into two beams with equal intensity at 90° polarisation angle difference. A quarter-wave plate is placed after the beam recombination to generate circular polarised light. Interferograms are taken at different transmission angles of a linear polariser just before the CCD/CMOS camera (Supplementary Fig. 8a). Here, we used a three-bucket temporal phase-shifting equation. Measuring the phase difference for samples of known ΔC, the relationship between ψ and ΔC can be established to obtain the contrast factor CF, which relates refractive index variation to concentration variation. CF is dependent on a number of factors including temperature, concentration, species and laser wavelength. Under isothermal condition for a single-wavelength PSI, same ΔC in different species has different contribution to ψ and this is why single-wavelength PSI cannot be used to determine the concentration of different species in the solution. Nonetheless, PSI is highly sensitive for detecting concentration differences in binary solutions. More details are available in Supplementary Method 3.

ICP measurement of ionic concentrations

Concentrations of cations in seawater substitute were measured in two laboratories at i) the Australian National University (ANU) using inductively coupled plasma mass spectrometry (ICP-MS, iCap RQ quadrupole, Thermo Fisher) and ii) ALS Limited, Environmental Division in Sydney, using inductively coupled plasma atomic emission spectroscopy (ICP-AES). At ANU, samples were gravimetrically diluted using 2% HNO3 with five replicates. The concentrations were drift-corrected using a gravimetrically measured external standard.

Molecular dynamics simulations

MD simulations were performed by Nanoscale Molecular Dynamics (NAMD) simulation software, following our previous work on non-equilibrium MD framework for assessing thermodiffusion in aqueous electrolytes52. A multi-component seawater brine solution was constructed as a model of natural seawater, containing the seven largest monoatomic ions that constitute seawater, that is Cl, Na+, Mg2+, Ca2+, K+, Br and Fl. The composition is detailed in Supplementary Table 1. More complex ions were not considered (i.e. SO42−) due to incompatibility with the simulated water. The water model used is TIP3P-FB. A control NaCl brine solution was constructed with the same molar concentration as the seawater brine, that is approximately 14,000 water molecules and 650 ions. The simulation volume was a rectangular prism of 20 × 5 × 5 nm with periodic boundary conditions, simulating a bulk solution. A quasi-linear temperature profile was maintained across the system by two spatially defined thermostats, spanning a temperature range of 20–60 °C. More details are available in Supplementary Method 6.

Three independent replicates of the NaCl brine and four independent replicates of seawater brine were simulated. After simulation calibration procedures, the ions freely diffused across the temperature axis in a constant volume and energy ensemble for 280 ns. The initial 20 ns is attributed to the time taken for the system to reach quasi-steady state, based on the scaling law53, and is excluded from final calculations. The time-averaged steady state concentration of ions in each brine solution was measured across the temperature axis with a bin spacing of 1 nm. Since MD is a discrete system, the instantaneous concentration within the small bin volume can fluctuate significantly. Thus, the quasi-steady-state system should be sampled for a sufficient number of times to obtain a statistically meaningful time-averaged concentration. Therefore, the simulations continued to run after 20 ns until 280 ns to allow a long enough time for adequate sampling to be performed. In Supplementary Figs. 9 and 10, it is shown that 260 = 280 − 20 ns yields a sufficient sample size. Ion concentration drop between the hot and cold boundaries was measured from a linear fit of the concentration profile against average temperature and used to characterise a Soret coefficient for seawater substitute brine and NaCl brine.