Salinity-Dependent Interfacial Phenomena Towards Hydrovoltaic Device Optimization

Evaporation-driven fluid flow in porous or nanostructured materials has recently opened a new paradigm for renewable energy generation. Despite recent progress, major fundamental questions remain regarding the interfacial phenomena governing these so-called hydrovoltaic (HV) devices. Together with the lack of modelling tools, this limits the performance and application range of this emerging technology. By leveraging ordered arrays of Silicon nanopillars (NP) and developing a quantitative multiphysics model to study their HV response across a wide parameter space, this work reveals the complex interplay of surface-charge, liquid properties, and geometrical parameters, including previously unexplored electrokinetic interactions. Notably, we find that ion-concentration-dependent surface charge, together with ion mobility, dictates multiple local maxima in open circuit voltage, with optimal conditions deviating from conventional low-concentration expectations. Additionally, assessing the HV response up to molar concentrations, we provide unique evidence of ion adsorption and charge inversion for a number of monovalent cations. This effect interestingly enables the operation of HV devices even at such high concentrations. Finally, we highlight that, beyond electrokinetic parameters, geometrical asymmetries in the device structure generate an electrostatic potential that augments HV performance. Overall, our work, which lies in between single nanochannel studies and macro-scale porous system characterization, demonstrates that evaporation-driven HV devices can operate across a wide range of salinities, with optimal operating conditions being dictated by distinct interfacial phenomena. Thus it offers crucial insight and a design tool for enhancing the performance of evaporation-driven HV devices and enables their broader applicability across the salinity scale of natural and processed waters.


Introduction
Water evaporation is ubiquitous in nature, occurring spontaneously even without solar illumination.It enables continuous energy exchange through the water cycle and is a currently untapped renewable energy source [1][2][3][4] .As an example, the total power generation potential of natural evaporation from lakes and reservoirs in the contiguous United States was estimated at 325 gigawatts, > 69% of the US electric energy production rate in 2015 3 .Various evaporationdriven devices have been explored which could be broadly classified as tandem devices 5,6 , hybrid systems 7,8 , and self-powered generators [9][10][11][12][13] .In particular, in 2017, among the self-powered generators, it was demonstrated that evaporation-driven flow through a functionalized porous carbon film could reliably generate sustained voltages up to 1 V and 100 nA current at ambient conditions 11 .Following this, various approaches such as a decrease of fluidic impedance, electrical resistance 14 , and surface functionalization 15 of the material structure have been demonstrated to improve the electrical power from 100 nW to ∼ 10 µW from a cm-sized device.
More recently, the use of Silicon has been shown to improve the current density significantly via ion-electron coupling through coulombic interaction, resulting in the enhancement in power density by two orders of magnitude 12,[16][17][18] .Overall, these results have attracted growing interest in the field of evaporation-driven HV devices 1,12,13,[19][20][21] .However, due to the complex micro-and nanostructure of the devices investigated so far, there remains a lack of fundamental understanding related to the underlying interfacial phenomena.Furthermore, the concentrationdependent performance of HV systems remains largely unknown.Interestingly, though, a recent study has shown that exploiting the synergistic effect of a nanofluidic diode and redox reactions 22 a substantial amount of energy can be harvested in high salinity conditions.Thus, it would be extremely appealing to explore the possibility of designing HV devices that can be operated optimally across salinity scales of natural water (fresh and seawater) and processed water (brine).
On the other hand, pressure-driven flow in functionalized nanochannels has been widely studied for electrical energy conversion via the streaming current pathway [23][24][25][26] .In particular, ion transport in nanopores with a diameter spanning from a few microns to sub-nanometers has been carefully studied in the Iontronic community [27][28][29][30][31] .This has shown that surface charge, pore geometry, ionic mobility, etc., are critical for the device output.In contrast, almost all the previous reports on HV used bulk porous material 6,9,11,12 , which limits the understanding of the role of geometrical parameters of the nanostructured material on energy conversion.
Additionally, previous studies on HV devices focused on a narrow parameter space, which could not highlight the relevance of many interesting phenomena observed in other electrokinetic devices, such as regime of EDL overlap, concentration dependent surface charge 25,27 , non-linear electrostatic screening 32 , effect of ionic mobility, and charge inversion 33,34 .Therefore, there remains a need to understand the fundamental role of a variety of electrokinetic phenomena in evaporation-driven HV devices.In particular, the interplay of structural, interfacial, and fluidic properties must be clarified in order to identify opportunities and guidelines for future device design and performance optimization.
Here we present a controlled study of evaporation-driven HV devices that reveals distinct interfacial phenomena across salinity scales and quantifies their impact on performance.By using a nanostructuring approach applied to arrays of Silicon nanopillars (Si NPs), we systematically change the solid-liquid interfacial area, the liquid confinement size as well as the ion concentration and type.We then correlate their effect to the open circuit voltage (VOC) and power output (Pmax) of the device.Importantly, by developing a quantitative Multiphysics model to interpret the experimental results, we provide deeper insights into the dominant solid-liquid interfacial phenomena at different regimes.For the first time, we show that chemical equilibrium at the interface plays a critical role in all tested conditions.In fact, it controls the dissociation of surface groups, leading to an electrolyte concentration dependent surface charge 35 .Additionally, we highlight that intrinsic geometrical asymmetries in the device structure boost the voltage output independently from electrokinetic modulation.In terms of concentration-dependent performance, we confirm that, at low concentrations (< 1 mM), geometry must be optimized to exploit electrical double layer overlap.Most interestingly, at intermediate concentrations (< 0.1 M), we demonstrate that, due to the interplay of concentration dependent surface-charge and ion mobility, multiple VOC local maxima exist, allowing for optimal operating conditions that deviate from low-concentration expectations.Additionally, by assessing the HV device response at high concentrations (up to 4M), we uniquely show evidence of ion adsorption and charge inversion for several monovalent cations (Li + and K + ).We also report viable power levels under these extreme operating conditions.Overall, our controlled nanostructuring approach, which lies in between nanofluidic studies and macro-scale porous system characterization, highlights the importance of controlling electrokinetic and geometrical parameters in HV systems.Our results also show that, by leveraging distinct interfacial phenomena, evaporation-driven HV devices can be developed for a broader range of concentrations, and hence applications, than previously thought.Therefore, in combination with the presented predictive model, this work opens new avenues for the engineering of evaporation-driven HV devices with improved performance and expanded applicability across the entire salinity scale of natural and processed water, from freshwater to brine.

Results and discussion
Understanding the fundamental mechanisms of voltage and current generation in evaporation-driven HV devices requires control over the solid-liquid interface properties and the liquid nanoconfinement, i.e., nanochannel geometry.Our devices consist of cm-scale regular arrays of Si NPs etched in a p-type Si wafer (Figure 1A).Using a combination of colloidal lithography and metal-assisted chemical etching (see Methods), we fix the pitch (p = 600 nm) of the hexagonal array of NPs (see red line in Figure 1C) while varying their diameter (Dnp) and length (Lnp) in the range 420 nm -560 nm and 1.23 μm -4.4 μm respectively.The corresponding mean pore diameter Dp ranges between 140 nm and 280 nm (Figure 1D and Figure S1).Therefore, changing the Si NPs dimensions directly controls the nanochannel geometry and the solid-liquid surface area (AS-L), as needed.
For hydrovoltaic testing, the sample is placed inside an HV cell and then wetted with 150 μl of deionized water containing different salts (KCl, LiCl, CsCl) of varying concentrations (from 1 μM to 4 M) while being placed in ambient conditions (T = 22-24 o C; and humidity = 25-30%).Thus, evaporation occurs naturally (Figure S2).The electrical response is measured using an Ag/AgCl electrode, placed in the liquid right above the Si NPs, and an aluminum contact, previously deposited on the back surface of the Si wafer (Figure 1A and Methods).During the electrical measurements, the Ag/AgCl electrode and silicon substrate do not contact each other, and the electrical circuit is complete as soon as the liquid is dispensed, thereafter voltage and current can be measured.Figure 1B shows a zoom-in view of a single nanochannel within our device.Upon wetting, the native silicon oxide layer on the surface of the Si NPs 36 will dissociate, resulting in a net negative surface charge σ (yellow circles), for pH above the isoelectric point. 35Concurrently, positive ions will adsorb on the surface (red circles) and an electrical double layer (EDL) will develop within the liquid.The distribution of ions in the confined liquid is governed by the Nernst-Planck-Poisson-Boltzmann equations (see Methods, eqns.(M1)-(M3)).The steady-state concentration distribution of the ionic species depends on different modes of ion transport, namely (i) conduction along an electrostatic potential gradient, (ii) diffusion due to a concentration gradient, and (iii) advection induced by the evaporation-driven flows (convection).
Any resulting imbalance in the ion distribution along the nanochannel leads to the measured electrical potential difference (Figure 1F).Overall, the total electrochemical potential of the i-th ionic component along the nanochannel length depends on the electrical potential, the ion concentration, and the chemical potential (see Methods eqn.(M4)).Interestingly, we observe that, due to the closed bottom surface of the nanochannel (z = 0), the studied geometry presents an intrinsic asymmetry in the surface-charge distribution.As a result, even in the absence of an evaporation-induced flow, an electrostatic potential, and therefore a VOC can be measured between the top and bottom surfaces of the nanochannel.Non-zero convection changes the distribution of ions and, depending on the flow profile, it can enhance or suppress the VOC.
Overall, our device structure presents three main advantages: (i) control over the liquid nanoconfinement, (ii) geometrical asymmetry within the nanochannel, and (iii) easy electrical contacting.Upon wetting, the VOC evolves with time starting from V(t 0 ), which is zero (negative value) for moderate (high) concentrations and attains a stable steady state value V(t ∞ ) (Figure S3).The VOC as well as the current-voltage curves (I-V) of the device (Methods and Figure S3) were recorded as a function of the physical and chemical properties of the system and will be discussed extensively later.

Modelling Approaches
The development of a suitable model for these evaporation-driven HV devices is essential to interpret the experimental results and unravel the complex electrokinetic interactions.Based on the aforementioned observations, the studied system can be described with the equivalent electrical circuit shown in Figure 1E.It consists of three parts: 1) a current source, representing the evaporation-driven streaming current Istr, and the related ionic resistance, Rionic, 2) a capacitor with a resistor in parallel, representing the EDL capacitance, CEDL, and the associated charge transfer resistance, Rct, due to geometrical asymmetrical and 3) an external load resistance, RL.
One can express the measured voltage at a given load resistance as follows: For the open circuit condition (RL → ∞), we thus obtain: Concurrently, the power to an external load can be calculated as   =   2 /  and the maximum power output (dP L )/(dR L ) = 0 is thus: The streaming current and the ionic resistance can be calculated from the ion concentration and the flow velocity (see Methods, eqns.(M5) and (M6)).However, no analytical expression is known for EDL capacitance for the geometry shown in Figure 1B.Therefore, we also developed a 3D numerical model (COMSOL®) to solve the Nernst-Planck-Poisson equation to determine the equilibrium distribution of ions and resulting electrostatic potentials (see Methods).To perform these calculations, however, we first need to identify an equivalent simplified geometry.
Considering a top-view of the Si NPs array (Figure 1c), we observe that it is possible to define a hydraulic diameter (Dh) for the resulting nanochannels, with Dmin as the minimum separation between adjacent NPs as well as AS-L, solid-liquid interface area normalized by sample area as: The 3D array of nanochannels can thus be reduced to an equivalent annular nanochannel geometry of inner and outer diameters R1 and R2, respectively, such that Dh and AS-L are equal to those of the original array (Figure 1C, D, and S1).A spatially uniform surface charge is used as a boundary condition.Importantly, as described by the Grahame Equation, its magnitude is governed by the diffuse layer potential and therefore by the chemical equilibrium between the silica surface and the electrolyte 35 (see Methods or Figure S4 ).This means that σ depends on the EDL thickness and therefore the bulk electrolyte concentration.The evaporative flux was instead used to define the mass flow rate in the channel.The computed potential difference between the top and bottom of the nanochannel is equal to the VOC and can be directly compared to experimental results (Figure 1D, F).

Role of surface-charge and ion transport at low-to-moderate concentrations
We first measure VOC as a function of ion concentration by using KCl solutions ranging from 1 µM to 0.1 M (Figure 2A).Indeed, this parameter plays a critical role both on the physical extension of the EDL and the chemical equilibrium at the interface.Regardless of the geometrical properties of the sample, our measurements clearly show a non-monotonic VOC trend as a function of electrolyte concentration.In particular, we observe a VOC peak for intermediate concentration values (0.1 -1 mM), which can be more or less pronounced depending on the sample.This result is in sharp contrast to the monotonic decrease in VOC predicted by the mean-field Poisson-Boltzmann theory with a constant surface charge condition 25 (Figure S5).Instead, using our COMSOL model that accounts for the chemical equilibrium at the solid/liquid interface, it is possible to correctly capture the experimental trend (Figure 2B).In particular, for a fixed value of Γ in the chemical equilibrium equation, we retrieve the appearance of the VOC peak in the intermediate range of concentrations.Interestingly, the model reveals that the VOC peak magnitude increases for increasing the length of the Si NPs and for decreasing Si NPs size (i.e., increasing nanochannel size), consistent with experimental results.Additionally, increasing the value of Γ from 8 to 12 results in a more pronounced VOC peak.This is due to the increase in the number of available surface sites that can readily dissociate to enhance the surface charge density as schematically shown in Figure 2C.Overall, the observed trend for VOC is explained based on the interplay between the magnitude of the surface charge density and the extent of electrostatic screening (Figure 2C), quantified by the Debye length (  ).This is defined as λ D 2 = ϵk B T /4πe 2 C b and is inversely proportional to the square root of the bulk ionic concentration (Cb).Thus, it monotonically decreases with increasing concentration.On the other hand, because of the chemical equilibrium at the interface, the magnitude of surface charge density and the diffuse layer potential both increase with concentration (Figure S4).Hence, a relatively high surface charge density and a moderate electrostatic screening give rise to an optimum condition and a peak in VOC at intermediate concentrations.To elucidate the significance of surface charge, we measured VOC for a single device at 1 µM concentration of KCl but varying the pH of the electrolyte across the isoelectric point by adding different concentrations of HCl (see time trace in Figure S6).As shown in Figure 2D, we observe the decrease in the magnitude of VOC with pH and then the change in sign below the isoelectric point.This highlights that the VOC is directly correlated to the surface charge, as the silanol groups at the surface could exist in one of these (Si-O -, Si-OH, SiOH2 + ) forms depending on the pH, thus modulating the surface charge density and sign.The values of Voc for deionized (DI) water are also shown in Figure 2A (left-most point in graph).DI water in standard condition has a hydronium ion concentration of 10 -7 M and therefore the Debye length is large compared to 10 -6 M KCl.Thus, from the perspective of EDL thickness, we might expect that DI water should give the highest VOC.Interestingly, we observe that, in the case of DI water, VOC is much less than for 1µM electrolyte.This observation is consistent with previous reports although an accurate explanation has not been provided yet.From our equivalent circuit model, we note that the open circuit voltage is directly proportional to the ionic resistance (Eqn.( 2)).Rionic, in turn, is inversely proportional to the mobility of ions.Thus, higher (lower) mobility ions are expected to result in lower (higher) open circuit voltages.In particular, for the studied system with a negative surface charge, we anticipate a stronger dependence on the cation mobility.Thus, in the case of DI water, the hydronium ions dominate the ionic transport, which results in much lower ionic resistance.These have much higher mobility due to the hopping mechanism 37 of ion transport compared to the potassium cations used in the electrolyte, leading to a sharp decrease in VOC.
To gain a deeper understanding of the effect of ionic mobility, we selected three chloride salts, namely LiCl, KCl, and CsCl, whose cations have increasing mobilities (Li + < K + < Cs + ) of 3.6, 7.32, 8.2 x10 -8 m 2 /Vs in aqueous medium 23 .For each salt, we varied the concentration from 1µM to 0.1M and recorded the VOC using the same device (Figure 3A).Both our experimental data and numerical modelling results (Figure 3B), which are in excellent agreement, clearly show that VOC increases with decreasing ion mobility, for the studied range of concentrations.Interestingly, we observe that the ionic mobility also influences the VOC peak at intermediate concentrations, which we previously related to the concentration dependent surface charge.We want to further highlight an interesting observation that connects the dependence of ionic mobility and the length of NP.The increase in length of NP and decrease in ionic mobility, both result in more pronounced intermediate peaks, which increase with surface charge, as shown in Figure 2A, B and Figure 3A, B. In particular, the lowest mobility cation (Li + ) exhibits the strongest VOC peak.
Next, to quantify the difference in ion mobility as well as the contribution of cations and anions to Rionic, we employed electrochemical impedance spectroscopy measurements for 1mM concentration of LiCl, KCl, and CsCl as well as KCl, KBr, KI. Figure 3C shows the diffusivity, and therefore the ionic mobility, determined by fitting the Nyquist plots using an appropriate electrical circuit (Figure S7).Interestingly, the large variation in diffusivity for different cations compared to the case of different anions confirms that cations dominate the ion transport in the studied system, as expected for a system with negative surface charge.
Next, we measure the I-V characteristic of the device (Figure S8) and determine the peak power (Pmax) at each concentration (Figure 3D).We observe a non-monotonic power generation dependence on concentration, having a minimum of around 0.1 mM concentration.Based on the previous discussion, the increase in Pmax at higher (above 1 mM) concentrations is attributed to the increase in Rionic, which dominates over the decrease in VOC as expressed in Eqn.(3).As a reference, we measured the bulk electrolyte conductivity for 0.01 mM, 1 mM, and 10 mM, as 1.64 μS, 95.4 μS, and 0.89 mS respectively.Importantly, we note that Rionic has contributions from both the bulk and surface conduction, which leads to a non-monotonic concentration dependent change in net ionic conductivity.This leads to an increase in power generation at low concentrations, which is attributed to the EDL overlap condition, where Rionic decreases due to an increase in surface conductivity.The relative contribution of the bulk and surface conductivity is governed by the concentration and size of the liquid confinement, which will be discussed in more detail in the last section.

Effect of high concentration (>1M) on surface charge and ion transport
The Poisson-Boltzmann equation is based on mean field theory where one can ignore ionion correlations as well as the adsorption of ions to the surface 38,39 .It thus holds for low electrolyte concentrations (up to ~0.1 M), when l B 3 C b ≪ 1, where l B = e 2 /εk B T the Bjerrum length of the solution and Cb is the bulk ion concentration, and ε is the dielectric constant of the electrolyte.However, for higher concentrations, the effect of ion-ion correlation becomes important.As a result, contrary to the Debye theory discussed above, the screening length (  −1 ) increases with concentration and becomes dependent on the size of ions.In particular, for concentrated electrolytes, this dependence scales as 32 , where a is the diameter of the counter-ion.Furthermore, at concentrations above ~1M, the influence of ion adsorption becomes relevant, leading to a reduction in the effective surface charge.This effect therefore reverses the trend discussed earlier of increasing surface charge with increasing concentration.
In particular, in the case of multivalent ions (Z ≥ 2), screening by counterions not only reduces the effective surface charge but can also reverse its sign.
This counterintuitive phenomenon is called surface charge inversion and arises due to the counterions forming a strongly correlated liquid (SCL) 34,40,41 , where the electrostatic interactions are much larger than thermal fluctuations (kBT).Surface charge inversion has been already observed for Z = 3 and 4 34 .Using pressure-driven streaming current measurements, it was also reported for Z = 2.However, there was an order of magnitude deviation from SCL theory in terms of the predicted minimum required ion concentration 33 .This deviation can be attributed to the omission of screening effects, which significantly enhances charge inversion 42 .Recently, direct measurements of surface excess charge using fluorescein particles have shown that at very high concentrations (>1M) even monovalent salts (Z=1) in silica confinement can give rise to charge inversion 43 .We thus explored experimentally the HV response of our devices for high electrolyte concentrations, including chloride salts of different-sized monovalent cations (Li + , K + , and Cs + ).
We show in Figure 4A the time traces of VOC for a few representative samples tested with 1M concentration KCl (purple and yellow curve), LiCl (green curve), and CsCl (pink curve).Upon wetting, we observe that all the curves exhibit an initial negative   0 value, where   0 =   ( = 0), which was not observed for lower concentrations (see the gray curve for 1mM KCl concentration, and Figure S3).Interestingly,   0 becomes more negative with increasing salt concentration (from 0.5M up to 4M) as well as for decreasing the length of the Si NPs and increasing liquid nanoconfinement (Figure 4B.left).This can be directly related to the effect of counter-ion adsorption close to the Stern-plane leading to a reduction, or even inversion, of the surface charge.We note here that, based on Eqn.M8, a change in Voc sign must be related to surface charge inversion.Yet, a positive Voc cannot exclude such an effect due to the additional voltage term introduced by the asymmetrical device structures.
Considering a fixed surface charge of 1e/nm 2 , which was measured previously for silicon HV structure 12 , we can estimate the magnitude of the short-range interaction, which can give rise to charge inversion, as 2.7, 4.5, 8.8, and 14.3 kBT for ion valence Z=1, 2, 3, and 4 respectively 41 .
This suggests that for monovalent cations, the high-concentration regions near the surface are mobile and are not expected to sustain charge inversion when equilibrium is achieved.Indeed, for some samples, the measured VOC changes over time, eventually stabilizing on a positive   ∞ value even at very high concentrations (Figure 4A, purple and pink curves, and Figure 4B (right)).In these cases, we can quantify the effective voltage generated by each sample as the difference between the initial and final VOC values, i.e.    =   ∞ −   0 .As a function of concentration, we observe a small but positive increase in    (Figure S9), consistent with the increase in screening length at high concentration.Interestingly, we also repeatedly observed steady-state sign inversion of VOC, which we postulate is related to charge inversion at high concentrations.(Figure 4A, green and yellow curves).
In order to overcome sample-to-sample variability and verify this observation (i.e., Figure 4A, yellow and purple curves), we first naively apply the SCL theory to calculate the critical concentration C0 34 .For Z=1 and σ ~1e/nm 2 we obtain C0 ∼ 1M.Based on this, we measured   ∞ for 4 different samples (S1-S4) at 1M electrolyte concentration using different cations like Li + , K + , and Cs + (Figure 4C).We observed that negative   ∞ , and therefore charge inversion, is most likely in the case of small cations like Li + , while no flipping in   ∞ sign was observed for big cations like Cs + .This result is consistent with previous direct measurements of surface excess of fluorescein 43 , which was positive for LiCl and NaCl (implying charge inversion of the silica surface) and the corresponding decay length at high concentrations.The dependence of ion size effect on charge reversal can be indeed attributed to EDL potential decay length at high concentrations, which scales as k −1 ~lB C b a 3 , and the specific adsorption that leads to excess cationic surface charge.Overall, in terms of HV device performance for operation in saline conditions, like seawater or brine, these results suggest that specific electrokinetic interactions should be accounted for while engineering surface charge and device geometry.For example, based on the composition of the natural water 41 , it would be important to engineer the interface properties to limit the surface charge density such that the critical concentration, C0 is much higher than the concentration of ions in water.Furthermore, together with the result of the previous section, these measurements confirm that effective operation in high concentration solutions (i.e.brine) would also be possible.

Effect of geometry and liquid nanoconfinement
The impact of changes in surface charge and screening length on the device performance are strongly dependent on the size of liquid nanoconfinement, and therefore the geometry of the device.Although we have shown some selected results in the previous sections, we now analyze more in detail the effect of the mean nanopore diameter, Dp, and length, Lnp.Firstly, we observe that the effect of geometrical parameters on measured VOC and Pmax is complex and nonlinear (Figure 5).Indeed, changes in Dp and Lnp modify the interfacial area (AS-L).This, in turn, determines the net surface charge, the streaming current (Istr), and the associated ionic resistance (Rionic, see Eqns.M4-5).In particular, based on Eqns.4A-C, AS-L decreases with Dp and increases with Lnp, leading to an overall increase in Istr.Additionally, we observe that Rionic can be also expressed as: Rionic ≡ 1/ (1/RB + 1/RS), where RB = Lnp /(πD p 2 σB) is the bulk ionic resistance and σB is the bulk conductivity, while RS is the surface resistance caused by the presence of the EDL and its distinct ionic conductivity.At small Dp (or Dmin), (Eqn.4C), the EDL overlap increases, resulting in a higher space charge density and therefore a lower surface resistance contribution to Rionic, as well as a higher streaming current.
Overall, the effect of surface resistance is quantified by a non-dimensional Dukhin number, Du = σs/(σBDp), and it has been shown to scale as 44 Du~ σλ D 2 /eD p .This shows that the role of surface resistance becomes prominent in the case of small Dp as well as for low ionic strengths.Hence, with the increase in electrolyte concentration, the bulk ionic resistance decreases, while the surface resistance can increase significantly, and therefore the relative contribution of the two is primarily determined by the diameter of the nanopore and bulk ionic concentration.From Figure 5A, B, we identify three different regimes with respect to change in the size of the nanoconfinement (i.e.Dp) in the following order: 1) EDL overlapping regime, where, for larger Dp values, the space charge decreases leading to a lower VOC (blue shaded area).
2) Streaming current dominated regime, where an increase in DP leads to higher streaming current (red shaded area).3) Ionic resistance dominated regime, where a decrease in ionic resistance limits VOC (grey shaded area) which can level off or even decrease.The concurrent variation in Pmax is due to the interplay between the ionic resistance and VOC as shown in Eqn.
(4).It can be noted from our measurements that not all samples/electrolyte concentrations show the existence of all regimes.This is due to the limited number of samples with different geometrical parameters used for measurements and some intrinsic variability in the initial surface charge.
Finally, with an increase in the length of the NP, the ionic resistance increases, and therefore VOC is expected to increase too.However, for a very long etching time, the surface composition can change, towards a reduction of the available surface sites on the NP due to HF removal of the native oxide layer on Silicon.As a result, we observe a VOC saturation for large etching times, i.e., for longer Si NPs, as shown in Figure 5C.We would furthermore draw the readers' attention that other factors such as a change in the fluidic flow pattern and sticking of long NPs due to surface forces can also affect the electrokinetic response.However, the discussion related to those aspects is beyond the scope of this manuscript.Interestingly, if we consider the peak power perunit solid-liquid interface area (Figure 5D) we observe that it decreases with an increase in interface area.This means that an optimum value of peak power exists for a given pitch and diameter of the NPs.Overall, this shows that geometrical parameter (Dp, AS-L, and LNP) plays a critical role in controlling the performance metric of the HV device which has a complex interplay with the surface charge.In particular, for optimizing the device performance one has to carefully identify the regime of operation for a given set of geometrical parameters and variability of surface charge depending on the fabrication process.

Conclusion
To summarize, we found that ion-concentration-dependent surface charge, together with ion mobility, entails optimal HV operating conditions that deviate from conventional lowconcentration expectations.Additionally, we observed that HV devices can successfully operate at concentrations exceeding the critical value (C0), charge inversion affecting the sign and magnitude of the steady-state Voc.This leads us to conclude that to-date studies of evaporationdriven HV devices have been unnecessarily limited to deionized water or low ion concentrations (regime of EDL overlap for ~10-100nm liquid confinement size).Instead, our results show that there is ample margin for improving their performance and extending their scope within the wide range of salinity conditions available in natural and processed water.More broadly, our controlled experiments and unique quantitative modeling clarify the complex interplay of surface charge, liquid properties, and geometrical constraints in evaporation-driven HV systems.They indeed uncover how different combinations of solid and liquid physical properties define distinct limiting regimes of operation.In this regard, Figure 6A highlights four main controlling quantities we have identified, namely ion concentration (C), surface charge (σ), solid-liquid interface (Asolidliquid), and liquid nanoconfinement (Dp), the latter two depending in our devices on the length (Lnp) and diameter (Dnp) of the Si NPs.It also shows how changes in the magnitude of these parameters, impact Voc and Pmax.Importantly, within this complex and multidimensional parameter space, it is possible to take advantage of different interfacial phenomena and identify suitable operating conditions for different salinity values.We found that for fresh water conditions (~1 mM), EDL overlaps is necessary.This can be realized, as expected, under low total surface charge and small nanochannel size (high confinement).Yet, by increasing the surface area or the surface charge, larger nanochannel sizes become viable.Furthermore, under seawater conditions (from ~10mM up to ~100mM), thanks to the chemical equilibrium at the interface, by controlling the surface charge (i.e.Gamma value in Figure 2B) an optimum can be engineered at large Dp values (> 100nm).This implies that nm-scale geometrical confinement can be avoided, simplifying the scalability of HV devices.Finally, at high salinity levels (above 1M) charge inversion can be leveraged by minimizing the solid-liquid interface and initial surface charge.Yet, long-term operation at very high concentrations can be challenging due to ion adsorption and salt crystallization, which directly affect the surface properties and geometry of the nanostructure (Figure S10).Thus, it will require further investigation.
From the perspective of device geometry engineering, our results confirmed that structural asymmetries lead to an open circuit voltage contribution entirely due electrostatic effect.Importantly, this term is what generates a sizable Voc (>0.1 V) in our system, despite very low evaporation rates (Figure S2).This could be further enhanced by engineering a spatially nonhomogeneous surface charge distribution through chemical or physical processing 45,46 .
Interestingly, taking advantage of our comprehensive model for predicting the performance of HV devices, we expect that the VOC performance metric can be further augmented by improving the rate of evaporation.Figure 6B shows that, depending on surface charge and geometry of the nanoconfinement, Voc can be doubled by a five-fold increase of the evaporation rate compared to ambient conditions, which means that power can be increased up to four times.This is largely due to the enhancement in streaming current and to be confirmed it will require a more in-depth understanding of the fluid dynamics in HV devices.Overall, based on these guidelines, by knowing the detailed composition of ions in the used water, it becomes possible to optimize geometrical and interfacial properties of the evaporation-driven HV devices.Thus, our work, which lies between nanofluidic studies of individual nanochannels and macro-scale porous device testing, offers critical insight into how to enhance the performance of evaporation-driven HV devices and points towards broader application opportunities for these self-powered systems.

Fabrication of Silicon nanopillars array
Metal-assisted chemical etching (MACE) of crystalline silicon combined with colloidal lithography was used for the fabrication of a cm-scale array of silicon NPs 47,48 (Figure S11).It involves the selfassembly of polystyrene nanospheres at the water-air interface.Then the non-closed-packed assembly of PS nanospheres was compressed to a pressure of approximately 25-30(N/m 2 ) using the Langmuir-Blodgett system, which resulted in a homogeneous closed-packed hexagonal lattice of PS nanospheres 49 .The closed-packed monolayer was then transferred to Pirhannacleaned Silicon substrates which were diced into 2cm X 2cm chips.Plasma etching was used to reduce the diameter of the PS nanospheres with an initial diameter of d = 600nm.After goldsputtering deposition of thickness 20nm and lift-off, a gold nanomesh is formed which is used as an etching mask for MACE.Prior to gold sputtering, 3-5 nm of Ti is sputtered as an adhesion layer.
This forms a stable contact between gold and the substrate to avoid delamination during MACE.
The liftoff was done by putting the substrate in Toluene and ultrasonicating it at moderate power for 3-5 minutes at room temperature.Finally, MACE was performed by putting the substrate in an aqueous HF/H2O2 solution with a volumetric percentage of HF and H2O2 as 10% and 2% respectively.The diameter of the NPs is controlled by changing the time of plasma trimming of the polystyrene nanosphere, while the length of the NPs is controlled by changing the MACE time.The obtained sample was then treated with oxygen plasma (60 sec, 1000 W), to improve the hydrophilicity and surface charge.

Electrical Measurements
The electrical measurements configuration is shown in Figure 1A, in which silicon NP is the active substrate and the Al layer acts as a back electrode while Ag/AgCl is used as a top electrode.For comparison, we tested the planar silicon device and NP device with Ag/AgCl electrodes and different porous top electrodes (Figure S12).The planar silicon gives almost zero VOC.We used an Ag/AgCl electrode as it is considered fully reversible, which ensures that the charges accumulated in the electrode EDL are entirely consumed by the electrodes at overpotentials, which are practically zero, and there is no unwanted potential difference, which induces a conduction current.The open circuit I-V and EIS measurements were done using CHI bipotentiostat.During EIS measurement, for covering the entire area of the active Silicon NPs, a porous graphite electrode was used instead of Ag/AgCl.I-V characteristics were obtained using linear sweep voltammetry (-0.5-0.5V) at a scan rate of 0.01V/s.The EIS measurements were done at zero applied DC voltage with an amplitude of 10mV in the frequency range of 100Hz-1MHz.

Modelling and simulation
∇ ⋅   +  ⋅ ∇  = 0 () Performing a full 3D simulation for a hexagonal lattice is computationally expensive, so we simplified our simulation, by transforming it to an equivalent annular cylindrical geometry.
During this transformation, the two important parameters, the hydraulic diameter, and the solidliquid interfacial area were kept constrained which gives us unique inner and outer radii R1 and R2 respectively.An evaporative flux is used to impose the flow rate condition (Figure S2).A spatially uniform surface charge was used as a boundary condition.
The magnitude of surface charge was governed by the equilibrium between the silica surface and the electrolyte M6.In the case of an isolated surface with a curvature or non-overlapping double layer, the charge density on the surface satisfies the Grahame equation 35 M7, and by solving this equation coupled with the M6, we can obtain surface charge as a function of electrolyte concentration.In our model, we use Gamma as a free parameter as it can vary largely depending on the surface preparation.So here, we use different values of Gamma based on the available literature to show how the number of surface sites available for deprotonation affects surface change, and therefore the VOC.We did not attempt to retrieve the experimental curve using our simulation, because we used the Grahame equation together with chemical equilibrium boundary conditions that hold for isolated surfaces.No analytical solution is available for a nonisolated surface and in that case, one has to sweep sigma for a range of values of surface potential and then find the intersection point with the chemical equilibrium condition.
To give a theoretical description for geometrical asymmetry-induced potential difference we define the total electrochemical potential of the i-th ionic component, Φi, along the nanochannel length (z-direction), which is given by: where Ψ  is the electrical potential, determined by the overall charge distribution, where Φ 0 i , c 0 i , zi , µi, ex are the standard chemical potential, standard concentration, valence, and excess chemical potential 50   The electrolyte here is KCl in DI water.

S4: Surface charge modelling
The equations were solved as a function of concentration at pH=7, pK=7.2, and C=0.29 Fm -2 , ε=78.5, T=300 K, and average radius of curvature, r=300 nm.The inset shows the equivalent electrical circuit used to fit the results and obtain the required circuit elements.The W2 is the Warburg impedance, which can be used to calculate the diffusivities.The diffusion coefficient can be obtained from the following equation 2 : Figure S8: a) IV characteristics were obtained using Linear Sweep Voltammetry (LSV) by sweeping the voltage from -0.5V-0.5V.Using IV data, the maximum power output was determined.b) The IV characteristics for different concentration.Direct contact with electrochemical means that the Ag/AgCl is contacted directly with the silicon substrate, while in case of electrochemical contact, the Ag/AgCl electrode was put in the electrolyte 'film' on top of silicon substrate.In the Dry condition, the device behave as a diode with zero current at zero applied voltage.In presence of electrolyte, the curve shifts downwards having a negative short circuit current and positive open circuit voltage.For 1mM, we intentionally make a direct contact of Ag/AgCl wire with Si substrate, and we see a diode like behavior but with negative short current.The magnitude of current in this case is consistent with the increasing trend with concentration.The simulation was performed using Comsol 6.0 using the transport of dilute species, and electrostatics packages that were coupled.The simulation was performed in 2D-axis symmetric domain as shown in the figure above.The normal flux of ions at the wall was set to zero, while at the top a concentration boundary condition was used.A concentration dependent surface charge was used as the boundary conditions at the solid surface as described in S4, while potential was set equal to zero at the top surface.The computed potential difference between top and bottom is equal to VOC.The meshing was finer near the wall due to large gradients, while a slightly coarser mesh was used near the bulk to optimize the simulation time. is to walls is much higher than in the bulk, while it is opposite for anions.For high bulk concentration, the electrical potential decays much faster due to small Debye-length, and therefore the concentration of ions is more uniformly distributed (lower concentration gradient) throughout the channel.For low bulk concentration, the electrical potential decays much slowly away from the surface, and therefore the ions concentrations have larger concentration gradients.

Figure 1 :
Figure 1: Mechanism of the HV device made of Silicon NPs array.A) Schematic of the electrical measurement setup for Silicon NPs hydrovoltaic device using a single Ag/AgCl electrode.B) The cross-section of the nanoconfinement formed due to the separation between adjacent NPs.The surface contains a negative surface charge (yellow circles) that forms an EDL.The evaporation-driven flow causes the streaming of ions, Istr, and the inherent asymmetry causes a difference in EDL potential between the top and bottom surface even when Istr = 0. C) (Top) SEM image showing a top view of the fabricated array of silicon NPs.The red line shows the hexagonal arrangement of the NPs, while the blue line represents the triangular unit cell used for simulation.(Bottom) Cross-sectional view.D) (Left) 3D schematic of the triangular unit cell of the hexagonally arranged NPs showing the geometrical parameters of the nanostructures (pillar diameter, Dnp, pillar length, Lnp and mean pore diameter, Dp).(Right) Annular cylindrical nanopore geometry was used for simulations, including the calculated electrical potential distribution (see also panel (F)).E) Equivalent electrical circuit in which streaming of ions and the asymmetrical EDL potential between top and bottom are represented by a current source, and a capacitor in parallel with a resistor respectively.The ionic resistance and external load are represented by the appropriate resistances.F) Vertical cut-plane of the simulated cylindrical nanopore in panel (D) showing the counter-ion concentration distribution, with ion flux (left), and the electrical potential distribution, with electric field lines (right).The bulk ionic concentration is 10 µm KCl for the given simulation results.

Figure 2 :
Figure 2: Effect of Low to Moderate Electrolyte Concentrations.A) Dependence of the measured VOC on the electrolyte concentration for a series of samples with different diameters (Dnp) and lengths (Lnp) of the NP.Lines are only a guide to the eye.B) Dependence of the simulated VOC on the electrolyte concentration.Two different values of the free parameter gamma (the number of available sites on silica surface) are used to show the importance of surface charge modulation (solid lines Γ = 8, dashed lines Γ = 12).Lines are only a guide to the eye.C) Schematic of the individual NP, showing a qualitative difference in surface charge density for Γ = 8 and Γ = 12.The middle graph shows a qualitative difference in the surface potential ψ and decay of EDL potential due to a change in surface charge density.D) Measured VOC at 1μM concentration for different pH values.The flipping of the sign of the surface charge across the isoelectric point is reflected in the sign of the measured VOC.Depending on the pH, the magnitude of I -, I + , and I 0 changes that modulate the sign and magnitude of surface charge.

Figure 3 :
Figure 3: Effect of ionic mobility on HV device performance at low-to-moderate concentration.A) Measured concentration-varying VOC and its dependence on the counter-ion mobility.The studied electrolytes are LiCl (black curve), KCl (blue curve), and CsCl (green curve) in DI water.B) Calculated VOC as a function of electrolyte concentration for different ionic mobilities, following the order Li + < K + < Cs + .Two different values of the free parameter Γ are used to show the variation with respect to the surface charge density.C) Diffusivity/mobility (related by Nernst-Einstein relation) of salts with different anions and cations determined from electrochemical impedance spectroscopy measurements for a constant concentration of 1mM.D) Measured maximum power output as a function of electrolyte concentration (KCl) for different diameters and lengths of the NPs.The colors correspond to the same devices as in Figures 2A and 4B.

Figure 4 :
Figure 4: Effect of electrolyte concentration on the measured VOC at high concentration.A) Time-trace of VOC for various concentrations and geometrical parameters of the NP.In contrast to low concentration, the voltage at t=0 is negative instead of zero.B) Measured VOC at t=0 (  0 , left) and steady-state VOC, (  ∞ , right) as a function of KCl concentration and for different samples (lines are a guide to the eye).C) Measured   ∞ , for 4 different samples (S1-S4) of 1.2 μm length of the NPs and varying diameter of the NPs, at 1M concentration for different monovalent cations.Charge inversion (i.e.negative   ∞) is more likely to happen for small ions like Li + than Cs + .

Figure 5 :
Figure 5: Dependence of geometrical parameters on the measured electrical outputs.A) Measured VOC as a function of the mean pore diameter (Dp) obtained by varying the diameters of the Si NP.Data are for a fixed length of the NP equal to 2.3 µm.B) Measured peak power as a function of the mean pore diameter for the same length of NP and various concentrations.The different colored regions represent the different regimes of ion transport as described in the text.C) Measured VOC as a function of the normalized solid-liquid interface area, AS-L.This is varied by changing the length of the NP (Lnp), as shown on the top axis, for a fixed NP diameter equal to 480 nm.D) Measured power normalized by solid-liquid interface area.The respective lengths of NPs are shown on the top axis.

Figure 6 :
Figure 6: Strategies towards boosting HV device performance.A) Contour plots of the measured VOC and Pmax across a wide range of salinity conditions available in natural and processed water.At low salinity, optimal performance can be achieved by leveraging the EDL overlap.At intermediate salinity, optimal operation can be achieved by engineering the surface charge and geometrical parameters, which is also highlighted by the intermediate peaks in Figure 2A-B.At very high salinity, one can leverage the charge inversion regime for optimal operations of the device.B) Simulated open circuit voltage for a series of different geometrical parameters as a function of evaporative flux for two values of surface charge.This shows that the performance can be further augmented by fluid dynamic consideration and enhancing the rate of evaporation from the system.
)For modeling the evaporation-driven electrokinetic conversion phenomenon in an array of vertical NP, we considered different modes of ion transport using the Nernst-Planck equation for the transport of dilute species, coupled with the Poisson-Boltzmann equation for equilibrium distribution of ions.The details of the simulation such as package use, and boundary conditions are given in S13.Diffusion becomes relevant because only the bottom part of the nanopores has surface charge, streaming of ions due to evaporation-driven flow, and migration of ions due to EDL potential.Based on the potential and charge distribution, it is also possible to compute the streaming current and ionic resistance to be used in the equivalent electrical circuit model:

Figure S1 :Figure S2 :
Figure S1: The SEM image of different samples used in the experiments was used for analysis to obtain the mean pore diameter and porosity as a function of Nanopillar diameter.

Figure S4 :
Figure S4: Surface charge density and diffuse layer potential for two different values of the parameter gamma.The chemical equilibrium and Grahame equation were solved for a particular value of Gamma as a function of electrolyte concentration that gives a unique surface charge density and diffuse layer potential value. 1

Figure S5 :DFigure S6 :Figure S7 :
Figure S5:To elucidate the importance of concentration dependent surface charge, we performed simulation with a fixed surface charge of -10 mC/m 2 .This result, which are monotonically decreasing with concentration, is in sharp contrast with the experimental trends.The monotonic decrease is due to increase in surface charge screening characterized by the Debye length, while the surface charge density is fixed.

FigureFigure S10 :Figure S11 :Figure S13 :
Figure S9: a) Effective open circuit voltage.b) initial open circuit voltage for different devices and range of KCl concentration from 0.5-4M.

Figure S14 :
Figure S14: Cation and anion distribution for a negatively charged surface for a bulk electrolyte concentration of 1uM (left) and 100uM (right) KCl under purely static condition (with no convection).Concentration of cations close of the i-th component, respectively.At equilibrium, the total S1: SEM Image statistics