Water, not salt, causes most of the Seebeck effect of nonisothermal aqueous electrolytes

,

Industries discard thermal energy on a large scale, and tapping into this resource may help society with its muchneeded energy transition [1].Among the alternatives, electric and electrochemical devices with temperatureinduced open-circuit voltages are attractive as they have no moving parts [2][3][4].The generation of a thermovoltage ∆ψ = ψ hot − ψ cold by a device subject to a temperature difference ∆T = T hot − T cold is called the Seebeck effect, characterized by the Seebeck coefficient S = −∆ψ/∆T [3,5].The Seebeck effect of solid-state devices relies on electrons and holes moving apart in thermal gradients [6].Electrochemical cells filled with aqueous [7] and polymeric [8] electrolytes and ionic liquid/organic solvent mixtures [9,10] show a Seebeck effect as well, which is usually explained by anions and cations moving apart in thermal gradients.Faradaic processes can also cause a Seebeck effect in cells with redox-active electrolytes [5,11].
Ions move in nonisothermal fluids due to the interactions among themselves and with solvent molecules [4,12].On mesoscopic length scales, the thermodiffusion flux of an ion species i can be written as J th i = −Q * i ∇T , with Q * i being the ions' heat of transport and T being temperature.The Q * i of aqueous electrolytes relate to the Gibbs free hydration energy G hyd i and hydration entropy S hyd i [13,14] by First, S hyd i is negative (and Q * i positive) for most ions, so when electrodes block fluxes, cations and anions pile up near the cold electrode.This so-called Soret effect of spatially varying salt concentration can be probed through its impact on the electrolyte's conductivity [15][16][17] and refractive index [18,19], giving experimental access to Q * i .Second, the size and valency of ions affect their hydration shell and, thus, their hydration entropy.According to Marcus's theory [20], smaller ions have larger hydration shells and more negative S hyd i .Third, Q * i tends to be larger for cations than for anions [21].If so, thermodiffusion of an electrolyte between electrodes at different temperatures leads to excess cations near the cold and anions near the hot electrode.This ionic charge separation causes a potential drop between the electrodes: the Seebeck effect.Combining J th i and the Nernst-Planck equation for ionic fluxes due to diffusion and electromigration, an expression can be derived for the steady-state Seebeck coefficient [7,18,22,23].For binary electrolytes, where e is the elementary charge.As Eq. ( 2) followed from an extended Nernst-Planck equation, it accounts for ionic thermodiffusion and mean-field electrostatic interactions.Finite ion sizes and ion-ion correlations, important for dense electrolytes and ionic liquids [24][25][26], are ignored, so Eq. ( 2) may not hold for such fluids.
Recent experiments on nonisothermal cells filled with aqueous [7] and polymeric [8] electrolytes found Seebeck coefficients above S = 1 mV K −1 , much higher than predicted by Eq. ( 2).Another experiment by Xu et al. [27] found potential differences between blocking electrodes held at different temperatures in separate electrolyte reservoirs, ruling out thermodiffusion as its cause in this case.Together these studies suggest that mechanisms beyond thermodiffusion contribute to the electrolyte Seebeck effect.Molecular dynamics (MD) simulations may help identify these mechanisms.One MD study determined G hyd i from simulations of isothermal systems at different temperatures [13].Others simulated nonisothermal electrolytes [13,18,28] and found that the thermal polarization of bulk water contributes to the Seebeck effect [29].None of these works accounted for electrodessimulating open, periodic systems instead.However, the electrolyte/electrode interface may well contribute to the Seebeck effect.MD simulations revealed a potential drop χ of about −500 mV over water/vapor interfaces; χ varied with temperature by up to −1 mV/K [30][31][32][33].If the surface potential at electrode/electrolyte interfaces also depends on temperature, holding two electrodes at different temperatures would produce different surface potentials and a potential difference between them.Here, we study nonisothermal thermoelectric cells with polarizable electrodes to determine how the electrode/electrolyte interface contributes to the electrolyte Seebeck effect.
We simulated pure water and aqueous CsF, KCl, NaCl, and LiI between two flat parallel blocking electrodes made from three graphene layers (see Fig. 1).All MD simulations were performed with the June 2022 version of the Large-scale Atomic/Molecular Massively Parallel Simulator (lammps) [34].We used the constant charge method of the electrode package [35] to keep each polarizable electrode overall uncharged while allowing the charge on them to fluctuate locally [36].We used different force fields for the SPC/E water [37], graphene [38], and halide and alkali metal ions [39].Nonbonded interactions had a cutoff of 1.2 nm, and the shake algorithm held bonds and angles of water molecules rigid [40].A PPPM k−space solver with a relative accuracy of 10 −6 computed the long-range interactions [41].The system was periodic in the electrode plane, and a correction for slab geometries was used [42].Initial simulations cells were around 2 nm × 2 nm × 16 nm.A force on one of the electrodes equivalent to a pressure of 1 bar adjusted the electrolyte's density.After letting the system equilibrate for a few nanoseconds, the distance 2L between the electrodes was between 11.4 to 12.2 nm, and the position of the electrodes was fixed.Thermostats set the temperature of the electrolyte by global velocity rescaling with Hamiltonian dynamics [43] in two regions immediately next to the electrodes, both of width ϵ = 1 nm.Upon applying a temperature difference ∆T , a linear temperature with a slope ∆T /[2(L − ϵ)] developed between the thermostats in less than 1 ns.All simulations were performed with a time step of 1 fs over 30 ns.For better statistics, simulations were repeated up to eight times with independent starting positions and flipped thermal gradients (thus up to 16 simulations).
We first discuss cells with 1440 water molecules and no ions.We use a Cartesian coordinate system with x and y lying in the plane of the electrodes of surface area A; the coordinate z runs from −L to L between the electrodes.The water's partial charges cause a spatially varying charge density ρ(z) and local potential ψ(z), according to the twice-integrated Poisson equation [29,30], where ε 0 is the vacuum permittivity, and z 0 is an arbitrary reference point left of the cell where we set ψ(z 0 ) = 0. We sample ρ(z) from the MD simulations by time-averaging the partial charges in bins spanning the xy plane and 1 pm wide in the z-direction.
Figure 2 shows ψ(z) for a "cold" and "hot" isothermal system at 293 K [Fig.2(a)] and 373 K [Fig.2(b)], respectively.In both panels, ψ(z) varies strongly near the electrodes but not in the cell's interior.To characterize ψ(z), we divide the cell into two interfacial regions L − δ < |z| < L of thickness δ and a bulk region where |z| < L − δ.We choose δ so that the bulk region is |z| < 4 nm.From linear fits to the potential in the bulk ψ fit (z) we determine the bulk potential drop, ∆ψ bulk = ψ fit (L − δ) − ψ fit (δ − L).Both Fig. 2(a) and (b) have ∆ψ bulk ≈ 0 mV, so water molecules are not polarized in the bulk.We define the surface potential drop χ at each interface as the difference between the electrode potential and the average potential at the edge of the bulk region, χ ± = ψ fit (±L ∓ δ) − ψ(±L).We observe similar surface potentials at both electrodes |χ − | ≈ |χ + |, with χ cold = −391 mV and χ hot = −489 mV.As a result, the potential difference ∆ψ = ψ(L) − ψ(−L) between electrodes, which can be partitioned as ∆ψ = χ − + ∆ψ bulk − χ + , is roughly zero.In general, the temperature dependence of χ is non-linear.Figure S1 in the Supplemental Material (SM) [44] shows that χ decreases roughly linearly between 293 and 373 K with (χ hot − χ cold )/80 K = −1.22mV/K, slightly higher than reported values for the SPC/E water/vapor interface [30,31].
Figure 2(c) shows ψ(z) in a cell with the electrodes held at 293 K and 373 K, respectively.The respective surface potentials are similar to those in panels Figure 2(a) and (b) and, therefore, depend on the local temperature: χ − ≈ χ cold and χ + ≈ χ hot .As the two surface potentials no longer cancel, we observe a nonzero ∆ψ and a corresponding Seebeck coefficient S = −1.25 mV/K [45].The nonisothermal cell also has a small bulk potential drop ∆ψ bulk ≈ 1 mV, in line with previous MD studies [29,30,46] [47].To separate surface and bulk contributions to S, we introduce a corresponding surface Seebeck coefficient, S surf = (χ + − χ − )/∆T , and bulk Seebeck coefficient, S bulk = −∆ψ bulk /∆T bulk , where ∆T bulk = ∆T (L−δ)/(L−ϵ) is the bulk temperature drop.As ∆T bulk < ∆T , these definitions yield S ̸ = S bulk +S surf .Table I lists all the resulting coefficients.
We show in Fig. 3(a) the near-electrode potential profiles of Fig. 2(c) shifted to zero at the electrodes (dotted lines) and their underlying charge density profiles (solid lines).The peaks and valleys of the charge density are inversely related to those of the mass density (see Fig. S2 in the SM [44]): where oxygen dominates, the charge density is negative and the mass density high; where hydrogen dominates, the charge density is positive and the mass density low. Figure 3 shows that the charge density ρ(z) varies less at a higher temperature, changing the surface potential χ.
To understand the temperature-dependent surface potential drop of water, we inserted multipole expansions of the charge density, ρ(z) = M(z)− d dz P z (z)+ d 2 d 2 z Q zz (z)+ . .., into Eq.(3) to determine how monopolar (M), dipolar (P z ), and quadrupolar (Q zz ) terms contribute to ψ(z) [29,30,32,46,[48][49][50][51][52][53][54][55].However, multipole expansions are not unique [56]: the terms M, P z , and Q zz contain an arbitrary reference point z ref m [see Eq. (S2) in the SM [44]].Figure S3 in the SM [44] shows P z and Q zz 's con- tributions to ψ(z) for five different z ref m (M = 0 for water).These data differ much so the different terms in a multipole expansion give little physical insight.Nevertheless, Q zz 's contribution to the potential difference between two points is proportional to the difference in Q zz (z) between those points [56].As Q zz = 0 in the electrodes, the quadrupole density does not contribute to ∆ψ.Consequently, in water-filled cells, only P z (and M, if ions are added) contribute to ∆ψ.
To unambiguously characterize the water's orientation, we consider charge-free layers that we define in terms of the cumulative charge of N a atoms with partial charge q i and position z i along the z-axis.
Note that neither Eqs. ( 4) and ( 5) contains any arbitrary reference points.Figure 3(c) shows P 1 , P 2 , and P 3 on the cold (blue) and hot (red) sides of the cell of Fig. 2(c).P l takes positive and negative values, and has larger absolute values at the cold side.The figure also shows the sum of the dipoles, three times larger near the cold side, and the dipole moment P center of the rest of the cell, which is much smaller than P 1 , P 2 , and P 3 .Summing all P l yields a nonzero overall dipole moment and associated potential drop ∆ψ = l P l /ε 0 .Hence, Fig. 3(c) shows that an increase in temperature leads to decreased water ordering and, in turn, a net dipole and potential difference.
We now turn to electrolytes with 1332 water molecules and 24 ion pairs of CsF, KCl, NaCl, or LiI between electrodes held at 293 K and 373 K.The ion's sizes are set by the σ ± LJ parameter of the Lennard-Jones force field; Table I 2 of Marcus [20] into Eqs.(1) and (2).
Table I shows us the following.First, our MD simulations yield S of the same size but opposite sign as experiments [7].This discrepancy may be caused by the different water structure near different electrodes; the experiments of Ref. [7] employed titanium, whereas our simulations concerned graphite electrodes.Second, Table I shows that |S| ≫ |S bulk | for all electrolytes, |S| ≈ |S surf | for most electrolytes, and that adding ions to water leads to a nonsignificant decrease in S for all electrolytes.These observations suggest that the electrolyte Seebeck effect is more caused by interfacial electrolyte ordering than thermodiffusion.Third, for all ion pairs except CsF, S surf values lie around S = −1.25 mV/K of pure water.This does not mean, however, that the electrolyte Seebeck effect is caused only by interfacial water structure.Adding ions to water changes χ − and χ + , but usually by the same amount, so the shifts in the surface potential drops cancel in S surf = (χ − − χ + )/∆T .Figure S4 in the SM [44] shows that the CsF anomaly is caused by a strong Cs + localization in the second water layer near the cold electrode, not present near the hot electrode.Fourth, S Agar and S Marcus correspond to bulk ionic thermodiffusion, but neither agrees with S bulk .Marcus's theory predicts different sized ions to have different S hyd i , so thermodiffusion should contribute strongly to the Seebeck coefficient of electrolytes with small and large ions [20].We indeed see that S bulk decreases monotonously with ∆σ LJ ; Marcus's theory predicts a monotonous increase instead.The discrepancies between S bulk and S Agar and S Marcus may be caused by temperature and concentration dependence of S bulk .Previous MD simulations found S bulk to depend strongly on temperature; for LiCl, S bulk changed sign around 305 K [29].Moreover, Agar's and Marcus's data concerned electrolytes at infinite dilution, while our simulations discussed so far were at 1 M.In the SM [44] (see Fig. S5), we consider different salt concentrations and find that S bulk generally increases with concentration.However, even at the smallest concentrations studied there (0.5 M), discrepancies between S bulk and S Agar and S Marcus persisted.Finally, the small values of S Agar and S Marcus further support our conclusion that the large S in our simulations cannot be explained by ionic thermodiffusion alone.
In summary, we presented the first MD simulations of nonisothermal electrolyte-filled cells with explicit electrodes.In our MD simulations, (i) the Seebeck effect of aqueous electrolytes is caused more by interfacial electrolyte structure than by ionic thermodiffusion; (ii) adding ions to pure water changes the interfacial electrolyte structure [57], but usually to the same degree near the cold and hot electrodes, leaving S surf unaffected; (iii) current models do not capture the thermodiffusion of dense electrolytes.Hence, further work is needed to fully understand the electrolyte Seebeck effect.
So far, many experiments concerned polyelectrolytes, ionic polymer gels, and ionic liquids [4].Experiments, theory, and simulations are more likely to meet for systems with simpler components, for instance, aqueous electrolytes, only a few of which have been studied [7].Next, it would be interesting to go beyond our idealized flat-electrodes setup.An electrode's shape, surface defects, and functional groups should affect the nearby electrolyte structure.Tailoring these properties could boost a system's Seebeck coefficient, making it more attractive for applications.
SUPPLEMENTARY MATERIAL to: Water, not salt, causes most of the Seebeck effect of nonisothermal aqueous electrolytes Ole Nickel, Ludwig J. V. Ahrens-Iwers, Robert H. Meißner, Mathijs Janssen

SURFACE POTENTIALS OF WATER
Figure S1 shows the surface potential χ for SPC/E water for several isothermal systems between 293 K to 373 K. Similar to Chapman and Bresme [30] and Sokhan and Tildesley [32], where the surface potential is discussed for different water models at the water/vapor interface, we find that χ decreases monotonously and the slope is reduced at higher temperatures.Hence, the Seebeck coefficient becomes larger when only smaller temperature differences at low temperatures are used.We fitted the data in Fig. S1

MASS DENSITY PROFILE
Figure S2 shows the mass density profile of water held between electrodes at 293 K (left) and 373 K (right).Density peaks near both electrodes are more pronounced near the cold than the hot side.Furthermore, the bulk density decreases from 0.993 g cm −3 to 0.950 g cm −3 in proportion to the local temperature, demonstrating the effectiveness of our density equilibration by applying a pressure of 1 bar to one of the electrodes to adjust the density.numerical issues).For unusual z ref m choices such as the hydrogen positions or the oxygen position mirrored across the hydrogen atoms, ψ P (z) has an opposite sign compared to the more common choices of the oxygen position or the centroid.The reason for these changes is that shifting the reference in the direction of the dipole will move molecules across the integration limit z in Eq. (S3a) depending on the value of their dipole.In fact, one can show that for shift ∆z ref m in the reference, a term proportional to ∆z ref m j q mj z mj is added to Eq. (S3a) and subtracted from Eq. (S3b).As there is no obvious choice for the reference point, ψ P (z) and ψ Q (z) carry little informative value.

IONIC CHARGE DENSITY IN NONISOTHERMAL CELLS
Figure S4 shows the charge densities of just the ions ρ ions (z) neglecting the water in a nonisothermal cell.All electrolytes have charge density maxima at different positions.The location of these maxima depends on the cation and anion's preferred location in the water layers, which is mostly based on the size of the ion and its hydration shell.Maxima signal that the cation is preferentially located there and minima mean that the anion is located there.Near the cold electrode, KCl, NaCl, and LiI are less layered than CsF.This explains their relatively small effect on the surface potential.CsF layers strongly at the cold side, with caesium ions located in the second water layer at the water's mass density peak (blue vertical line).The strong charge oscillations of CsF near the cold electrode cause a large increase in χ − compared to the surface pure water.CsF is much less layered near the hot interface, explaining the difference between χ − and χ + for this electrolyte.From the singular ions densities, we have seen that both cations and anions move towards the cold electrode, leading to a pileup in its vicinity as expected.

Figure 1 .
Figure 1.Simulation snapshot of 1 M LiI in water between graphene electrodes.Thermostats control the temperature in the blue and red regions independently.(Li + -yellow, I −purple, hydrogen -white, oxygen -red, carbon -cyan).

Figure 2 .
Figure 2. Potential profiles ψ(z) of water at 293 K (a), 373 K (b), and with a temperature gradient, where the left side is at 293 K, and the right side is at 373 K (c).

Figure 3 .
Figure 3.Further analysis of the water-filled cell of Fig. 2(c) with temperature varying from 293 to 373 K. (a) Charge density and local potential, where the hot side at 373 K is mirrored around z = 0, and its potential in vacuum shifted to 0 V.The vertical black line indicates the position of the first graphene layer.(b) Cumulative charge (lines) and oxygen number density (dotted lines).(c) Layer dipole moment of water in the first three neutral layers near the cold and hot sides of the cell, their sum, and the dipole moment of the bulk water in the remainder of the cell (black).

Figure 3 (
b) shows Q(z) for the nonisothermal cell of Fig.2(c).For this analysis, we moved the two hydrogen atoms of each water molecule to the point between them.In this way, Q(z) is a multiple of the oxygen charge and the points where Q(z) = 0 are easy to identify.For a single layer of water molecules, seen as a peak in the oxygen density, Q(z) rises, drops and passes zero, and rises again.Accordingly, we set the boundaries ζ l−1 and ζ l of uncharged layers l = 1, 2, 3, . . .at the second crossings of the cumulative charge, Q(ζ l−1 ) = Q(ζ l ) = 0, and the boundary of the first layer at ζ l−1 = z 0 .With this definition, the layer boundaries (vertical dotted lines) are close to minima in the oxygen number density (dotted lines).Figure3(b) clearly shows three layers near both electrodes; the boundaries of further layers are difficult to determine due to a weak signal-to-noise ratio.Accordingly, we only consider the first three layers near both electrodes and combine the rest of the cell in one big central layer.The distribution of atomic charges in layer l is captured by its layer dipole moment, FigureS1shows the surface potential χ for SPC/E water for several isothermal systems between 293 K to 373 K. Similar to Chapman and Bresme[30] and Sokhan and Tildesley[32], where the surface potential is discussed for different water models at the water/vapor interface, we find that χ decreases monotonously and the slope is reduced at higher temperatures.Hence, the Seebeck coefficient becomes larger when only smaller temperature differences at low temperatures are used.We fitted the data in Fig.S1by a quadratic function,

Figure S1 .
Figure S1.The surface potential of the SPC/E water/graphene interface at different temperatures.

Figure S3 .
Figure S3.Potential profiles for the same system as in Fig. 2c in the main text.Contributions to the potential from (a) the dipole moment, (b) the quadrupole moment, and (c) their sum.The contributions are shown for five choices of the reference point z ref m : the position of each atom z ref m = z ′ mj with j ∈ {0, 1, 2} and with the absolute atom positions z ′ mj = z ref m +zmj, the centroid of the three atom positions z ref m = j z ′ mj /3 and the oxygen mirrored across the hydrogen atoms z ref m = −z ′ m,oxygen + j∈{1,2} z ′ mj .These data correspond to a 10 ns simulation, where the position of every atom is sampled every 1 ps to calculate Pz(z) and Qzz(z).

Figure S5 .
Figure S5.Plot of S bulk for different electrolytes at different concentrations.

Table I .
Potential drop and Seebeck coefficients of water and several electrolytes in different cell regions (electrodes held at 293 K and 373 K, respectively).Mean values and uncertainties are determined by block averaging 16 simulations.
[20]ar et al.[21]reported Q * i and Marcus[20]S hyd i for OH − and H 3 O + , but we did not include these ions in our MD simulations.
lists the differences in ion sizes ∆σ LJ = σ − LJ − σ + LJ for each electrolyte.The table also presents MD results for the potential drops χ − , ∆ψ bulk , and χ + , which we use to determine the Seebeck coefficients S, S surf , and S bulk .Last, Table I contains predictions S Agar and S Marcus , for which we inserted Q * i data from Table 1 of Agar et al. [21] and S hyd i data from Table