Electrochemistry, ion adsorption and dynamics in the double layer: a study of NaCl(aq) on graphite

Graphite and related sp2 carbons are ubiquitous electrode materials with particular promise for use in e.g., energy storage and desalination devices, but very little is known about the properties of the carbon–electrolyte double layer at technologically relevant concentrations. Here, the (electrified) graphite–NaCl(aq) interface was examined using constant chemical potential molecular dynamics (CμMD) simulations; this approach avoids ion depletion (due to surface adsorption) and maintains a constant concentration, electroneutral bulk solution beyond the surface. Specific Na+ adsorption at the graphite basal surface causes charging of the interface in the absence of an applied potential. At moderate bulk concentrations, this leads to accumulation of counter-ions in a diffuse layer to balance the effective surface charge, consistent with established models of the electrical double layer. Beyond ∼0.6 M, however, a combination of over-screening and ion crowding in the double layer results in alternating compact layers of charge density perpendicular to the interface. The transition to this regime is marked by an increasing double layer size and anomalous negative shifts to the potential of zero charge with incremental changes to the bulk concentration. Our observations are supported by changes to the position of the differential capacitance minimum measured by electrochemical impedance spectroscopy, and are explained in terms of the screening behaviour and asymmetric ion adsorption. Furthermore, a striking level of agreement between the differential capacitance from solution evaluated in simulations and measured in experiments allows us to critically assess electrochemical capacitance measurements which have previously been considered to report simply on the density of states of the graphite material at the potential of zero charge. Our work shows that the solution side of the double layer provides the more dominant contribution to the overall measured capacitance. Finally, ion crowding at the highest concentrations (beyond ∼5 M) leads to the formation of liquid-like NaCl clusters confined to highly non-ideal regions of the double layer, where ion diffusion is up to five times slower than in the bulk. The implications of changes to the speciation of ions on reactive events in the double layer are discussed.


Introduction
][10][11][12] Graphite provides a model substrate to study such interfaces; current technological applications involving carbon materials e.g., activated carbon, which offers a high surface-tovolume ratio for maximum ion adsorption, are oen made up of more-or-less ordered microdomains of graphite. 13Improving the design of these materials, therefore, requires an understanding of the graphite-electrolyte interfacial structure and dynamics, where many details remain unresolved. 14t is at the interface between solids and liquid solutions where changes occur to the average structure and dynamics of charge carriers (cf. in the bulk) when the two phases come into contact, even in the absence of a potential bias.This results in the formation of the so-called double layer and an electric potential difference (drop) in this region.Early models of the double layer sought to explain the capacity for charge storage at the interface, which can be measured as the differential capacitance, C d .
6][17] At low electrolyte concentrations and with low applied surface potentials, this model effectively predicts C d , in simulations 18 and in systems containing metals, 19 even if the simplied interfacial geometry it implies is unphysical for the electrode/electrolyte interface. 20he value of C d measured at graphite is much lower than expected from consideration of GCS theory alone.3][24] Assuming a GCS electrolyte structure, C d was then used to estimate the electronic DoS of graphite. 23,24t moderate to high electrolyte concentrations-relevant to many technical applications-experimental observations suggest that the GCS model inadequately describes the solution-side structure of the double layer at graphite.6][27] This deviation is consistent with an observed inection point in the screening length of NaCl(aq) solution as a function of concentration around 0.5 M from surface force balance experiments. 28These changes have been attributed to ion-specic adsorption, 25 with possible partial charge transfer; 29 however, without a detailed microscopic description of the interface, the origin of the variations in solution behaviour remain unclear.
1][32][33] Nonetheless, the complexity of the system favours the application of atomistic simulations to fully characterise the double layer structure that results from changes to bulk electrolyte concentrations and applied surface potentials.These methods allow for a molecular description of how the complex interplay between ion adsorption, solvation and speciation in the double layer, coupled with any changing mobilities for charge carriers, affects the electrochemical and thermodynamic properties of the interface.
Calculations at the level of density functional theory (DFT) with implicit solvents have provided signicant insight into the nature of charge screening in the double layer. 29,346][37][38][39][40][41][42] Based on DFT calculations, classical pairwise interaction potentials were parameterised; [37][38][39] these capture the polarisability of the solution and carbon (in the form of ion-p interactions) at the interface, which can play a signicant role in structuring the double layer. 43,44Such classes of models further indicate asymmetric adsorption of ions, with cations likely to adsorb in preference to anions in the case of graphite.This contrasts with metal-electrolyte interfaces, where ion-p interactions are absent and anion adsorption likely dominates. 45The implications of these effects on the structure and dynamics of the double layer at graphite will be addressed in the present article by combining our own simulations and experimental measurements.
Constant chemical potential MD (CmMD) 46 simulations were performed to explore the double layer in the (electried) graphite-NaCl(aq) system.CmMD mimics open boundary conditions; thus, maintaining a constant thermodynamic driving force for ion adsorption at graphite, and conserving electroneutral solutions beyond the double layer.With this approach, we are able to relate the spatial extent of the double layer to the nature of charge screening in this region.We explain changes to the electrochemical properties of the double layer by relating the screening of charge to the potential difference across this region.In addition, we obtain a detailed description of the local (electro)chemical potential, speciation and mobility of ions orthogonal to the surface.The results indicate a signicant departure from ideal solution behaviour in regions conned to the double layer even at moderate levels of NaCl(aq) concentration in the bulk.Our work provides much needed molecular-level insight into the structure and dynamics of electrolyte solutions in contact with carbon surfaces over a wide range of concentrations.

Results and discussion
2.1 The structure of NaCl(aq) solutions at the graphite surface Simulation cells were prepared where a graphite slab, comprising eight graphene layers, was positioned at the centre of the simulation cell x axis in contact with NaCl(aq), such that the system was symmetrical about x ¼ 0. All methodological details are provided in ESI Section A. † CmMD simulations were performed using the GROMACS (v 2018.6)MD package 47 with the Plumed (v 2.5) Plugin 48 for a range of xed bulk solution concentrations: c b NaCl ¼ 0.23-1.05M and 1.2-9.2M in two system set-ups.The upper limit here signicantly exceeds the solubility of NaCl in water, and would not be possible to prepare experimentally, but is instructive to study computationally.The nanosecond timescales associated with MD give rise to vanishingly low probabilities for crystal nucleation and allow the metastable solution state to be investigated.Within minimal uctuations, the CmMD method successfully maintained the bulk concentration of cations and anions beyond the interface (see ESI Section C †).
Concentration proles orthogonal to the surface.Concentrations as a function of x for species in the solution phase are provided in Fig. 1 and also in full in Fig. S5 and S6.† Preferential adsorption of Na + was observed at the graphite surface over the entire concentration range sampled, in line with the predicted order of ion adsorption energies using the adopted force eld. 37olvated Na + ions directly coordinate to graphite, shown by the narrow peak at x ¼ 1.5 nm in Fig. 1A.Maximum concentrations exceed those in the bulk by approximately two orders of magnitude even at the lowest c b NaCl (0.23 M).An adjacent layer of Cl À ions was observed, separated by 0.1 nm at the peak concentrations, and in line with other simulations. 49The Cl À concentration in this layer also signicantly exceeds the bulk value, although its broader width highlights a diffuse ordering of the anion.Considering the adsorbed layer of Na + to represent an effective surface charge density, this picture of charge screening is qualitatively consistent with the predictions of the GCS model; with the Cl À diffuse region of the double layer representing the counter-ion charge in solution.However, no obvious boundary between this region and a diffuse layer is apparent in the concentration proles, and any binary assignment of surface-bound states neglects the complexity of the dynamic adsorption in the rst ion layers.
As the bulk concentration of ions is raised (see Fig. S5 †), a clear departure from the above screening behaviour is observed.Around c b NaCl ¼ 0.5 M, the Cl À peak narrows and a second, diffuse layer of cations around x ¼ 2 nm emerges that exceeds bulk ion concentrations.Fig. 1A highlights that at c b NaCl ¼ 1.05 M, the concentration of cations exceeds that of anions at x ¼ 2-2.4 nm, and a hierarchical ordering of ions with opposing charge is apparent.A more compact double layer region is evident at the highest concentrations (see Fig. 1B and  C, with the complete data set available in Fig. S6 †).In the range c b NaCl ¼ 1-9 M, a shi in the position of the Cl À rst peak by Dx z À0.03 nm highlights a contraction of the rst ion layers, and the diffuse Cl À peak at 1.2 M was resolved into two clear peaks (with a second peak emerging at x ¼ 1.83 nm at around 6 M bulk ion concentrations).This is concomitant with a shiing of the second Na + peak away from the graphite surface, which ultimately splits into a rather diffuse doublet peak which connes the Cl À layer.
The crowded structure at higher concentrations is reminiscent of the double layer in molten LiCl at planar electrodes under the constraint of a constant applied potential. 50Indeed, in double-layer capacitors containing ionic liquids, steric crowding at the electrode is a common feature, 51,52 which has been conrmed by atomic force microscopy (AFM). 53,54The similar response of the double layer structure to changes in bulk concentrations and applied surface potentials was identi-ed in simulations over 40 years ago. 18As well as the increasingly non-monotonic concentration proles (shown clearly in Fig. S7A † on increasing c b NaCl ), a shi in the ratio of maximum Na + : Cl À concentrations occurs around 3 M due to the narrowing of the second Cl À peak (see Fig. S7B †) which represents a signicant departure from the double layer structure at the lowest concentrations.
Perturbations to the solvent structure were apparent, mainly due to the presence of graphite.Three peaks are observed (see Fig. 1A) in the water concentration proles at the lower end of the bulk ion concentration range simulated.Solvent layers are separated by 0.3 nm; this distance was determined for water at clean graphite surfaces in a recent study combining simulations with AFM measurements. 55Only limited ordering of the orientation of water molecules perpendicular to the graphite surface was observed.This is highlighted in Fig. S5, † which shows that the peaks for water O and H atoms are approximately at the same position in x.Maximum values of c Cl at the interface tended to be found where the densities of water oxygen atoms are close to a minimum.Additional ion layers at the highest concentrations induce additional complexity to the water structure in the double layer (see Fig. S6 †).
Extent of the double layer.The edge of the double layer region on approach to the surface was marked by the position where solutions deviate from electroneutrality; hence, hc Na (x)i s hc Cl (x)i (where angular brackets indicate the mean concentrations in 0.5 nm moving windows in x).The size of the double layer from this position to the minimum in x for the solution phase is provided in Fig. 2A; this indicates that the double layer contracts as the ionic strength in the bulk solution initially increases, reaching a minimum around (c b NaCl ) 1/2 ¼ 0.76 M 1/2 (c b NaCl z 0.6 M).Beyond this bulk concentration, the double layer size increases, plateauing around 2 nm above (c b NaCl ) 1/2 ¼ 2.4 M 1/2 , with some noise in the data.Overall, using this composition measure, the double layer size is 0.6-2.2nm.The bulk electrolyte concentration marking a transition in the double layer structure (0.6 M) and the size of the double layer, within the concentration range sampled here, are in very close Fig. 1 Ion and water molar concentration (c) profiles from CmMD simulations as a function of x: the distance from the centre of the simulation cell.The graphite carbon surface (with zero applied charge) is shown by the grey peak on the left of the x axis.(A) Provides Na + (blue), Cl À (red) and water oxygen (O wat ; black) concentrations when the concentration of ions in the bulk is 0.23 M (solid line) and 1.05 M (dashed line).(B and C) Provide the cation (top panel) and anion (bottom panel) concentrations as a function of x in simulations targeting the higher end of the total concentration range.The colour scale here indicates the mean molarity of ions in the bulk.The shaded dashed lines mark the maximum in the first sodium peak (grey) and first two chloride peaks (green) at the highest bulk concentration.Arrows are provided to highlight the changes in the profiles as ion concentrations in the bulk are increased.A 0.3 nm excluded region separates the edge of the graphite basal plane from the ionic solution; this is due to atom centres being used to calculate the concentration profiles.
agreement with the results from surface force balance measurements which determine the screening length in NaCl(aq) solutions between parallel plates. 28he Debye length (k À1 ) is the characteristic length over which the electrostatic effect of a charge carrier in solution decays, and is proportional to 1/(c b ) 1/2 .This is derived from a linearised Poisson-Boltzmann equation (see Section B † for details), and is assumed to accurately determine the size of the double layer at low electrolyte concentrations.Fig. 2A indicates that k À1 decreases monotonically as c b NaCl increases.At relatively low concentrations, one can reconcile the decreasing double layer size from simulations by considering that increasing charge densities close to the surface will lead to a less diffuse double layer as contraction of the layers occurs.At high concentrations, however, the ionic crowding near the surface induces further perturbations to the solution away from the interface, and the double layer size increases with c b NaCl .A theoretical framework to predict the 'capacitive compactness' of the double layer was recently presented; 56 this indicates the dependency of the size of the interfacial region upon the ion valency.
It is instructive to consider the size of the interface region where the solvent structure is perturbed (cf. the bulk), which turns out to be independent of c b NaCl over the entire concentration range sampled.ESI Section D † details these measurements which indicate an interface region that is 1.4 AE 0.3 nm in size, that is approximately 3-5 water layers from the graphite surface.The structuring of water at planar interfaces appears to be rather insensitive to the electrolyte concentration and, for most practical purposes, the substrate material and surface contamination. 55ean ion concentrations in the double layer.NaCl concentrations at the interface (c i NaCl ) can be measured by integrating (c Na (x)c Cl (x)) 1/2 in regions of the proles in Fig. 1.To ensure a fair comparison between different cases of bulk concentration, NaCl , the concentrations of ions at the interface are three times greater than those in the bulk and the decay in the relative interfacial concentrations is proportional to (c b NaCl ) À0.25 (see Fig. 2B).While nite ion size effects clearly play a role in the local ion concentrations in the double layer, the total concentrations of ions in this region vary continuously with bulk concentration and can therefore be predicted without the need for simulations at specic concentrations.It is important to note that the total concentrations of cations and anions over the entire double layer region are equal, as shown in Fig. S8, † and that signicant ordering at the highest bulk concentrations means that, locally, ion concentrations can signicantly exceed the bulk (see Fig. S7A †).
Electrical properties of the double layer.Even at uncharged graphite surfaces, asymmetries in the adsorption of ions leads to deviations from local electroneutrality, as shown by uctuations in the solution charge density, r, as a function of x.Following the Poisson equation, this leads to varying electric elds, E, and electric potential, j orthogonal to the surface.In the above equation, 3 is the permittivity of the medium (3 ¼ 3 0 3 r , where 3 0 and 3 r refer to the permittivity of vacuum and the relative permittivity, respectively).Fig. S9 † provides these quantities over the full concentration range.These indicate that, at the limit of large x, E(x) and j(x) converge to zero, corresponding to the solution bulk.
A screening factor, f, can be dened as, with n indicating ion number densities; n Na appears in the denominator as cations consistently adsorb in the rst solution layer next to the substrate.Fig. 3A provides f(x) for the entire concentration range sampled.At the lowest c b NaCl , f(x) increases smoothly and converges to one in the solution bulk, consistent with screening by a diffuse anion layer.As c b NaCl increases beyond around 0.5 M, and the compensating anion charge layer becomes more compact, over-screening of the cation charge occurs and f(x) > 1. Over-screening in molten salts is a phenomenon that has been known for some time. 57In ionic liquids at electried interfaces, over-screening was suggested as a possible control on the electrochemical kinetics at the interface. 58he electric potential in x due to the charge distribution of ions (j i ) can be calculated using eqn (1) and (2), noting that r i (x) ¼ e(n Na (x) À n Cl (x)), where e is the elementary charge: The relative permittivity of the solution medium as a function of x (3 t r (x)) is affected by the proximity of interfaces [59][60][61][62] and ion concentrations. 63A recent simulation study from Olivieri et al. 64 indicated that truncation of the long range water molecule dipole moment correlations due to the presence of the substrate and the electric eld anisotropy at the interface leads to increasing 3 t r away from the surface.We adopted this model for the near-linear increase of 3 t r , from a value of around 10 in the rst solvent layer up to 71 (the bulk value of 3 r for the extended simple point charge water model 65 in our simulations) over x z 4 nm, as described in ESI † Section E-where we also discuss alternative models and their implications for the calculated j i (x) values.
The j i (x) curves calculated using eqn (3) are provided in Fig. 3B.When c b NaCl T 0.5 M, the sign of j i alternates due to the crowding of ions in the double layer.The value of j i at the graphite surface (with zero applied charge) is called the potential of zero charge (PZC): j PZC i .This potential difference is shown in Fig. 3C as a function of the bulk concentration.It is clear that the effect of increasing c b NaCl is to decrease j PZC i with a À0.020 AE 0.003 V M À1 gradient at moderate bulk concentrations.An inection point is observed when c b NaCl z 6 M, where further increases in bulk ion concentrations result in positive changes to j PZC i .The slope, (dj PZC i /dc b NaCl ) T,s is related to the so-called Esin-Markov coefficient. 66This is one of the few conventional means to experimentally assess the extent of specic adsorption, and would seem to provide an ideal way to relate the simulations to experiments measuring C d .Oen, C d data are recorded at widely spaced potentials, limiting the accuracy with which the minimum in C d -oen taken to represent PZC-is known. 25,27e therefore measured C d of freshly exfoliated HOPG as a function of potential over a range of concentrations, with 10 mV potential resolution, and examined the minimum in C d as a proxy for the PZC (j xpt ; see Section A2 †).With this ne potential resolution, the C d -j xpt curve displays two minima within 300 mV; the global minimum becomes deeper and shis to more negative potentials with increasing concentration (Fig. S2 †).Note that this feature remains visible at 50 mV potential resolution (Fig. S3 †), but would be lost at lower resolutions, oen reported in the literature. 22,23,25,27omputational models considering the effect of asymmetric ion adsorption indicate a shi to the PZC on incremental changes to bulk solute concentration that follows the sign of the preferentially adsorbing ion. 50However, these oen do not consider the situation of electrolyte solutions in which both cations and anions have favourable, but varying strength of interactions with the surface.Experimental studies sometimes attribute C d values solely to the adsorption of ions whose sign is opposite to the sign of the potential change relative to the C d minimum. 25This is a rather simple interpretation under conditions where the bulk concentration is far beyond the levels where alternating layers of charge emerge in the double later structure.The change in the value of j xpt identied at the C d minimum on increasing c b NaCl (see Section A2 † and Fig. 3C) is in good agreement with the mean change in j PZC i (c b ) found from our simulations, with a moderate negative gradient (À0.02V M À1 ) at lower concentrations.This contrasts with the linear shi as a function of ln (c b NaCl ) seen on Hg electrodes. 67The agreement provides considerable support to the simulation results.
The implication of eqn ( 3) is that increasing levels of ion screening, with the same underlying cation density distribution, results in positive changes to j PZC i .Conversely, if the screening of charges in the double layer is unchanging and the bulk concentration increases, then j PZC i becomes more negative.In our simulations, Fig. S11 † shows that at the limit of large x, NaCl , the resulting response to the PZC can inform about the structure of the double layer.

Ion activities in the double layer
The excess density of ions in the double layer has implications for the chemistry of this region.It is helpful to evaluate, therefore, contributions to the electrochemical potential, mĩ, of ions, i, as a function of x: here, m 0 is a reference chemical potential; a i (x) and z i are the activity and valency of species i, respectively; and, b ¼ 1/k B T, where k B and T are the Boltzmann constant and temperature, respectively.The total m̃for NaCl(aq) in x can be written as, where u is the fraction of cations and a AE (x) is a positiondependent average activity.We dene a AE (x) ¼ a Na (x) u(x) a Cl (x) 1Àu(x) , which converges to the mean ion activity, , in the electroneutral bulk solution where u(x) ¼ 0.5. 68Eqn (4) and ( 5) therefore reduce to the standard equations for the chemical potential (m) in a bulk homogeneous solution where j(x) ¼ 0 in the absence of an external electric potential.At equilibrium, a AE (x) and j(x) are stationary, and the electrochemical potentials across the double layer and in the extended solution are equal: Based on ts to the chemical potentials of ions explicitly calculated in their own simulations and from other in silico studies, 69,70 Zimmerman et al. 71,72 provided an analytical model to calculate the mean chemical potential of solvated ions using the adopted force eld: where b is the mean molality in units of mol kg À1 , and g AE is the mean activity coefficient for ions; therefore, a AE ¼ b NaCl g AE , and, where A ¼ 0.568 m À1/2 , B ¼ 1.17769 m À1/2 and C ¼ 0.177157 m À1 (where m ¼ mol kg À1 ).4][75][76] Values of ln g AE and mÑ aCl calculated using the above model are provided for the range of b NaCl in Fig. S12.† From CmMD simulations, the NaCl molality as a function of x is calculated from atom density proles according to (n Na (x) n Cl (x)) 0.5 /(0.018n wat (x)).Molalities in the bulk (b b NaCl (x)) were calculated from averages in the molality proles in stable regions far from the interface.b b NaCl (x) were substituted into eqn ( 8) and (7) to calculate (m NaCl À m 0 NaCl )/2 which equals the right hand side of eqn (6).Note that this approach assumes equal contribution of cations and anions to the mean ion electrochemical potential.The CmMD simulation technique used here ensures accurate estimates of the bulk chemical potential of ions-within the adopted model in eqn ( 7)-where cation and anion concentrations must be uniformly equal within a small uncertainty.
Analyses were performed using the density proles in Fig. 1, with Fig. 4A providing the NaCl molality proles for ve systems.These reach a maximum at x ¼ 1.6-1.7 nm: the position close to the minimum following the rst peak in cation concentration proles.An increasing peak around x ¼ 2 nm matches with the increase in cation concentrations in this region at high concentrations.Counter-intuitively, the region in x around the maximum b b NaCl (x) corresponds to a minimum in b À1 ln a AE (x) (see the shaded region Fig. 4B).This is due to a maximum in (2u(x) À 1)ej(x) as shown in Fig. 4C.The fact that at x ¼ 1.63 nm, u and j are at a minimum, results in the large positive contribution to m̃from this term in eqn (6).Beyond around 1 nm from the graphite surface, we nd that the contribution of (2u(x) À 1)ej(x) to m̃is zero, and in this region mr educes to m.
An important implication of the above observation is that the local ion molality (or concentration) at interfaces is not a good proxy for the electrochemical potential of ions.In systems where interfaces and extended liquid phases are in equilibrium, this is usually not problematic, due to the equality m̃¼ m b across the boundary layer and into the bulk solution.Although any partial charge transfer of surface-bound ions should be considered. 29,77In some surface-driven processes at equilibrium, knowledge of the interfacial structure might still be essential to predict outcomes.For example, in processes like salt precipitation, the rates for nucleation are affected by the kinetic factors associated with the supply of ions to growing crystalline embryos. 12Increased ion molalities close to the surface (being around ve times the levels of the bulk, when the bulk molality equals the equilibrium saturation level of 3.7 mol kg À1 ) likely mitigates the barriers to these processes.In non-equilibrium processes, knowledge of both the local molality of solute species and the electric potential is essential to determine m.

Ion correlations and diffusion
Simulations at the atomic level are perfectly suited to provide details regarding the collective arrangement and motion of ions in the liquid near the graphite surface.In this section, we characterise the species which emerge in the double layer at uncharged graphite surfaces and their diffusion.
Ion speciation.We investigate ion speciation by calculating the average rst shell coordination number, N i-j , between atoms i and j described in ESI Section A. † The average coordination numbers, N Na-Cl ¼ 0.02 AE 0.02; N Na-Ow ¼ 5.89 AE 0.05 and N Cl-Ow ¼ 7.23 AE 0.07 were calculated for a 1 M NaCl(aq) bulk solution simulated for 10 ns.The structure of the solvated ions agrees well with other simulation and experimental studies. 78A majority of ions form solvent shared and solvent separated ionpairs in the bulk (represented by the peaks at r z 0.5 and 0.7 nm in Fig. S13 †).
Fig. S14 † provides the average N as a function of x evaluated using Gaussian kernel (with 0.03 nm bandwidth) probability densities.These indicate that the average N Na-Ow (x) decreases by a value of one from the lowest to highest sampled concentrations in the bulk regions as N Na-Cl (x) changes from zero to one, with these changes becoming signicant when c b NaCl T 1 M. Interestingly, N Cl-Ow (x) is far less sensitive to changes in c b NaCl , which remain around the value identied in the bulk at 1 M over then entire bulk concentration range within the model.The N Na-Cl (x) coordination proles in Fig. S14 † A indicate a greater number of directly coordinated ions beyond the position of maximum densities in x for the rst ion layers.Essentially, the increased coordination occurs in regions of the double layer where there is a high density of both cations and anions.The maximum in N Na-Cl (x) shis to smaller values of x and an additional peak emerges at x $ 2.2 nm as crowding in the double layer increases.The features of these proles roughly correspond to the proles of b b NaCl (x) in Fig. 4A, which is a good indication that the association of ions in the double layer is due to increased ion densities.

At the highest levels of c b
NaCl , N Na-Cl > 2 at the maximum positioned at x ¼ 1.75 nm, as shown in Fig. 5A.Concomitant changes occur to N Na-Ow (x) in this region, with an additional minimum at x ¼ 1.6 nm.The exceedingly high anion concentrations here affect the ability for water to fully solvate cations.The increased Na-Cl coordination is due to the formation of many contact ion pairs that dynamically (dis)associate on the timescales of the simulations, leading to extended liquid-like networks of the type identied in the inset of Fig. 5A.These structures begin to emerge at the graphite surface in signicant numbers when the concentration of ions in the bulk exceeded approximately 5 M (i.e., beyond the nominal equilibrium saturation level for this force eld).The networks are reminiscent of other liquid-like ionic networks identied in simulation studies, 79 which were suggested as precursors to crystalline phases. 80,81o analyse these structures further, we performed cluster analyses using the method of Tribello et al. 82 (see ESI Section A † for details).Fig. 5B indicates that the effect of increasing c b NaCl (x) is to increase the number of ion clusters (dened as species containing more than two ions in direct contact) in the bulk and within the double layer.In the bulk, these clusters contain a maximum of three ions.In the double layer, however, the clusters contain more than ten ions at the highest c b NaCl (x), with a wide distribution in the size of the largest clusters due to the rapid time evolution of ion-ion correlations.Both the size and geometry of the networks rapidly changed over several nanoseconds of simulation and exchange of ions with the surrounding solution occurred.Similar liquid-like NaCl clusters were identied in simulations beyond the limit of solution stability (15 mol kg À1 ) in the bulk, where a change in the mechanism for salt precipitation occurs. 83Some experimental studies posit the existence of NaCl clusters even at moderate saturation levels. 84Further studies are now needed to explore the role that liquid-like clusters play on the nucleation of NaCl at interfaces.
Diffusion in solution.The diffusion coefficients, D, for ions were measured using the Einstein relation, described in detail in ESI Section A. † For reference, the average D for ions measured in simulations of bulk of 1 M NaCl(aq) was 1.14 AE 0.05 Â 10 À5 cm 2 s À1 .D and D x as a function of x are provided in Fig. S15; † these indicate that the surface decreases the diffusion of ions and water molecules within the double layer.The decrease in both D(x) and D x (x) on approach to the substrate is monotonic; hence, the diffusion coefficients closest to the graphite surface and in the bulk region have been plotted as a function of c b NaCl in Fig. 5C for simulations sampling the higher end of the entire concentration range.The values of D were found to decay following an approximately exponential trend: D ¼ D 0 e Àlcb , where l is the so-called decay constant.In the bulk, D 0 were 1.344 and 2.396 Â 10 À5 cm 2 s À1 and l were À0.236 and À0.386 for ions and water, respectively.A more negative l for water indicates that increasing ion concentrations retards the mobility of the solvent molecules moreso than solute ions.At the interface, however, l for ions was À0.335 (D 0 ¼ 0.819 Â 10 À5 cm 2 s À1 ): more negative than for water (l ¼ À0.221; D 0 ¼ 2.248 Â 10 À5 cm 2 s À1 ), which is most likely due to the increased concentration of ions in this region and the changes to the speciation of ions, discussed above.Around 1 M, D i /D b ¼ 0.64 and 0.87 for ions and water, respectively (where D i and D b are the diffusion coefficients close to the graphite surface and in the bulk).At the highest concentrations sampled, D i /D b ¼ 0.24 and 0.18 (for ions and water, respectively).The arresting of particle mobilities is largely due to decreased diffusion perpendicular to the interface.The top panel in Fig. 5 provides the x component of D for water and ions at the interface and in the bulk.At 1 M, D i x /D b x ¼ 0.35 for ions and this reduces to 0.03 at the maximum bulk concentration.This reects the high charge densities close to the graphite surface.In contrast, D i x for water molecules is unchanged compared to D b x at the lowest concentration, but beyond 9 M, D i x /D b x ¼ 0.09, and the high salinity interface retards the mobility of water molecules in x nearly as signicantly as for ions.
Ion transport properties in the double layer are oen assumed to match with those in the bulk, e.g., when calculating z-potentials using electrokinetic ow apparatus.Diffusion coefficients for ions in solution near the graphite surface on the order 1 Â 10 À7 to 1 Â 10 À5 cm 2 s À1 indicate an increased viscosity in the double layer caused by the changing solution densities in this region.While direct coordination of cations to the graphite was evident, no specic surface-site binding was identied, and diffusion was particularly limited perpendicular to the graphite surface.This picture is arguably consistent with the idea of a 'dynamic Stern layer'. 85,86However, this term is unhelpful, 20 failing to recognise the dynamic equilibrium between ions in the rst and adjacent solution layers.No clear boundary (slipping plane) between the diffusion of ions in a specically adsorbed layer at the surface and in the diffuse region can be identied from D or D x in Fig. S15.† We refer the reader to a recent monograph by Döpke and Hartkamp, 86 where these effects are discussed in the context of electrokinetic phenomena.

Graphite with applied electric charge
To consider the effect of applied electric elds, charges were uniformly distributed to the outermost carbon atoms at the graphite basal plane, discussed in detail in ESI Section A. † Negatively and positively charged surfaces were simulated with charge densities, s, in the range jsj ¼ 0.004-1 e nm 2 (see ESI Table S2 †).Our approach neglects the electronic response of the graphite to charging the material, which must be considered when comparing to experiments.This is reasonable, considering that DFT calculations indicate that applying an electric eld to graphite induces equal and opposite net excess s, centered close to the edges of a graphite slab, perpendicular to the eld direction. 21he concentration of ions beyond the double layer was maintained in CmMD simulations (see Fig. S16 †) where c b NaCl ¼ 1.05 AE 0.03 M. Assuming Poisson-Boltzmann behaviour, Grahame's equation (see eqn (13) in ESI Section B †) provides a direct relationship between (the effective) s and the potential change across the double layer, Dj.Even within the monotonic regime, Poisson-Boltzmann approximations can fail to accurately predict s. [86][87][88] These models are, therefore, unhelpful, particularly to determine the interfacial properties of NaCl(aq) on graphite at technologically-relevant electrolyte concentrations, as discussed below.
Structural asymmetries in the double layer.Ion number densities as a function of x in Fig. S17 † show that when increasingly positive surface charges are applied, Cl À ions were pulled closer to the graphite surface in a less diffuse anion layer.Concomitantly, the maximum density in the rst cation layer decreased and was shied further away from the surface.Despite this, a more compact double layer emerges due to increased cation densities in a second cation layer, shown most clearly at high values of s in Fig. 6A.When the surface charge was made incrementally negative, the rst two cation layers both increased in density and peaks in n Na (x) appear sharper.The maximum in the anion density decreased and was shied away from the graphite, although the anion layer remained diffuse, with the density of the tails in the distribution actually increasing (see Fig. 6B and S17 †).
Structural changes were observed in the solvent layers at the interface; these were particularly signicant when s was large, in line with simulations elsewhere studying electried planar interfaces. 89Fig. S18 † shows that when s was large and positive, water oxygen atoms were pulled closer to the surface and peaks in n Hw (x) appear sharper.More interestingly, when s was more negative than approximately À0.5 e nm 2 , a restructuring of water molecules was apparent at the interface, with a splitting of the rst peak in n Hw (x), as solvent molecules arrange their hydrogen atoms towards and away from the graphite surface.
Changes to the excess ion number densities (n i Àn ref i , where n ref i are the densities when s ¼ 0) as a function of s are reported in Fig. S19.† These highlight that the double layer undergoes an asymmetric enrichment of charge-balancing counter-ions and depletion of co-ions.The asymmetries can be quantitatively evaluated by considering the excess counter-ion concentration at the charged graphite surface.This was done by taking the difference in ion concentrations, Dc ¼ c A À c B (where A is the counter-ion type and B is co-ion), from the maximum in concentration proles where Dx < 1.5 nm in the MD data.Fig. 6C shows that Dc increases monotonically from a value of zero when s ¼ 0; this is qualitatively consistent with the predictions of the GCS model.At the positively charged surface, Dc is negative when s is small.The repulsion of cations in the rst ion layer (beyond the electrode) leads to a small decrease in the maximum anion concentrations in the rst anion layer and a small increase in the maximum concentrations in the second cation layer beyond the surface (see Fig. S19 †).These observations are beyond the predictions of simple mean-eld models and are due to the complex interplay of charge screening and volumetric constraints due to nite ion size effects at intermediate-high bulk concentrations.
Electrode charge screening.Asymmetries in the shape of the excess solution charge densities (dened as r 0 (x) ¼ r(x) À r ref (x), where r ref is the s ¼ 0 case, reported in Fig. S20 †) appear in the rst ion layer, with a doublet peak emerging at the negatively charged surface.Despite this, the features of the resulting E 0 (x) and j 0 (x) curves (calculated using the Poisson equation, and also provided in Fig. S20 †) are largely symmetrical when comparing the oppositely charged surfaces.Differences occur in the amplitudes of the uctuations in these curves.Fig. S22A † highlights a greater change to the electric potential across the double layer region, Dj 0 , in response to applying positive charges (cf.negative charges) at the graphite surface.The potential difference diverges when jsj ¼ 0.1 e nm 2 .
To further understand the above divergence, we calculated an electrode screening factor: Fig. S21 † shows that, considering only ions in solution, overscreening occurs at the positively charged surface due to the high density of anions in the rst solution layer at the highest s. f 0 then decays until it converges to a constant value far from the graphite surface.Small changes to s affect the maximum in f 0 (x) and, therefore, the gradient in the screening proles as f 0 converge to their bulk values.On the negatively charged surface, however, there is an initial over-screening that is compensated by charges in an adjacent solution layer, leading to a minimum in f 0 around 1.6 nm.We calculated Ð x 0 0 ð1 À f 0 ðxÞÞ (aer ensuring that all f 0 proles converge to a value of one).The resulting curves at the highest values of applied potential are provided in Fig. S21.† A positive (negative) shi in the converged value of Ð x 0 0 ð1 À f 0 ðxÞÞ is found at the positively (negatively) charged graphite surface.This is consistent with a more diffuse screening of the positively charged surface.When all solution atoms were included in the analyses, the same trends were observed, as highlighted by Fig. 6D; although the data here are noisy due to uctuations in the water density proles, particularly for small values of jsj.These changes correlate with the changes to jDj 0 j, presented in panel D of Fig. S22A.† The asymmetric accumulation of ions manifests in a greater capacity, C ¼ s/Dj 0 , for the surface to store ionic charge when negative potentials are applied.Fig. 6E shows that maximum and minimum C occur at small negative and positive values of Dj 0 .This is in very good agreement with the experimental C values (also show in Fig. 6E) calculated using the integral of C d (Dj 0 ) measured for the HOPG-NaCl(aq) system when c b NaCl ¼ 1 M (described in ESI A †).The asymptotic behaviour of C(Dj 0 ) is opposite to the trends at positive/negative potential observed in simulation studies of molten LiCl at atomically at electrodes; however, these electrodes are not related to any specic material. 50When Dj 0 was highly negative, C is up to 2 mF cm À2 greater than when Dj 0 was large and positive.
Despite the close agreement between the simulated and experimentally determined values of C, the differential capacitance data (C d ¼ ds/dDj 0 ; provided in Fig. 6F) show deviations around Dj 0 ¼ 0. Here, the mean simulated differential capacitance reaches a maximum around 11 mF cm À2 (shown inset of Fig. 6F, with the average over the full range of potential difference sampled in the experiments being 6.1 AE 1.2 mF cm À2 from simulations).3][24] Our results, therefore, indicate that both C d dl and C d q contribute to the measured C d at the chosen electrolyte concentration around Dj 0 ¼ 0. This was conrmed by evaluating C d q explicitly using 1/C d q ¼ 1/C d À 1/C d dl , as shown in Fig. S23, † which indicates that C d q is of the same magnitude as C d dl , on average, around Dj 0 ¼ 0. Increasing the potential difference results in rapidly increasing C d q , with a larger magnitude in the gradient (dC d q / dDj 0 ) on the anodic branch of the C d q (Dj 0 ) curve.This means that, at large values of s, C d dl dominates the measured C d at graphite.Indeed, the shape of C d q (Dj 0 ) matches well to experimental measurements of the 'quantum capacitance' at graphene electrodes, although the size of C d q is an order of magnitude greater here at similarly large values of applied potential. 90This is due to the proportional increase in the electrode (electron) density of states per unit surface area that occurs as the number of graphene layers are increased, as determined in recent simulations at the DFT level. 91The DFT calculations show that a C d q of around 10 mF cm À2 (at the C d q minimum) in 6-layer graphene is increased by an order of magnitude when jDj 0 j z 0.5 V, in close agreement to the evaluated C d q here.2][23][24] This was motivated by that fact that the minimum C d of graphite is approximately an order of magnitude smaller than the C d measured for metal electrodes at comparable conditions.It was, therefore, implicitly assumed that C d was dominated by C d q . 21,22,24Outside of the GCS model, however, there is no reason to assume that graphite electrodes accumulate solution charge in the same way as metal electrodes, such that they should have similar magnitudes of C d dl ; indeed, this was briey considered early on. 23,24ur estimates of C d dl are, to some degree, force eld dependent.At the concentrations adopted, however, it is likely that other force elds, which capture (at least implicitly) the polarisation at the graphite surface, will determine a similar value for C d dl due to the partial 'saturation' of the multi-layered double layer structure at high concentrations.In addition, the possible adsorption of airborne contaminants 27 at the surface, and the effect of steps and other surface defects in HOPG should be considered when attempting any comparison between simulations and experiments.Nonetheless, our results are consistent with a dependence of C d on the solution composition, reported here and elsewhere. 25,26Asymmetries in the tails of the C d curves are accentuated in the simulation results, and these are apparent also in experiments-highlighted more clearly by the changing gradients of the C d (Dj 0 ) curve in Fig. S22B.† This is a common feature in studies of ionic liquids, 58 where the anisotropic, bulky charge carriers saturate opposing surface charge densities differently.

Conclusions
Little is known about the effect of concentration on the structure and dynamics of electrolytes in molecular detail at carbon surfaces over a wide range of bulk concentrations and how these affect the electrochemical properties of the interface.Here, CmMD simulations were applied to simulate the graphite-NaCl(aq) interface in equilibrium with constant concentration, electroneutral bulk solutions sampling a range of electrolyte concentrations (0.2-9.2 M).Na + accumulation at the graphite surface in a dynamically adsorbed rst layer (as conrmed by analyses of ion diffusion coefficients) results in an effective positive surface charge density that is charge-compensated by the accumulation of Cl À ions.At the lowest bulk electrolyte concentrations, the concentration of anions decreases exponentially (exponential t R 2 ¼ 0.94 at 0.2 M), with the double layer thickness shrinking with increasing concentration.Above 0.6 M, however, a transition in the screening behaviour occurs, with alternating layers of cations and anions (up to four-ve solvent layers in extent) forming before the effective surface charge is neutralised, leading to an increasing double layer size with concentration.The concentration for this transition is in good agreement with experimental studies which found a minimum in the NaCl(aq) electrostatic screening length at 0.1-0.9M. 28 The crowding of ions, and increasing over-screening in the double layer region, manifests in a reduction to the change in the potential of zero charge with incremental changes to the bulk electrolyte concentration, which was conrmed experimentally.
Charging the graphite surface, when c b NaCl z 1 M, allowed for evaluation of the double layer differential capacitance (C d ).That the average over the same region of applied potential in simulations (6.1 AE 1.2 mF cm À2 ) and experiments (5.54 AE 0.60 mF cm À2 ) were of the same order of magnitude highlights a problem with previous analyses of the double layer capacitance (C d dl ) of the graphite-electrolyte interface, which assumed that a GCS model was applicable. 22,24We emphasise that parallels cannot be drawn between the classical picture of the compact layer capacitance at metal electrodes and the graphite electrode; this led us to reconsider the relative contributions to the total C d .Estimates of C d associated with the response of the graphite electronic structure to charging (C d q ), using our simulated values for C d dl and total measured C d , shows that C d q is of a similar magnitude to C d dl around Dj 0 ¼ 0 (around 10 mF cm À2 ), increasing substantially with the value of applied potential.This is in line with recent simulations of many-layer graphene. 91The signicance of this result is that C d (at moderate-to-high electrolyte concentrations) does not simply report on the density of states of the graphite material, as has become the accepted norm in much of the electrochemistry literature.
Evaluation of interfacial molalities showed increases of up to six times the levels of the bulk solution in regions conned to the double layer, where the solution mass densities were up to three times the levels of the bulk (1030 kg m À3 at 1 M).Because both cations and anions adsorb in the double layer, a complex solution structure emerges which affects the thermodynamic and aforementioned electrochemical properties of the interface.Despite the high molalities, the mean ion activities are reduced due to the large local electric elds present, emphasising the non-ideality of this region.Cluster analyses reveal the presence of enhanced clustering at concentrations above 5 M, which has implications for the rates of surface-driven reactions.The clusters are reminiscent of those identied in simulations at the limit of solution stability in bulk solutions ($15 mol kg À1 ), where phase separation becomes spontaneous. 83Given that the clusters emerge in the metastable solution (i.e., beyond the nominal solubility of NaCl in the bulk and below the limit of solution stability), further investigation is required to determine their role in salt precipitation facilitated by surfaces.
The insight gained here for planar carbon interfaces also applies to the wide range of activated carbon materials, which have graphitic microdomains.Further work will be required to understand how the structure at the planar interface is modi-ed by the complex, conned geometry and chemical heterogeneity of such materials.

a 1 .
5 nm solution region closest to the graphite surface was integrated, and the c i NaCl normalised by c b NaCl are provided in Fig. 2B.The plot shows a rapid decay in c i NaCl /c b NaCl on increasing c b NaCl , with concentrations at the interface converging to those in the bulk when c b NaCl > $6 M. At the lowest c b

Fig. 2 (
Fig. 2 (A) The depth of the double layer region, d, in CmMD simulations with varying scaled bulk ion concentrations, (c b ).The grey points provide the size of the double layer in x, with error bars indicating uncertainties of one standard deviation in the data.The black line provides a moving average for the data.Note that the minimum value of x for the solution phase (with reference to Fig. 1) was taken to be the mid-point between the carbon surface and the minimum value in ion concentration profiles: x z 1.35 nm.Also shown in the figure is the blue curve which is the Debye length (see ESI Section B †). (B) Provides the normalised, mean interface ion concentrations (c i NaCl ) in a 1.5 nm solution region close to the graphite surface vs. the bulk concentration of ions.The black line is a fit to the data with functional form 1.8(c b NaCl ) À0.25 and error bars show uncertainties of one standard deviation.

Fig. 3
Fig. 3 Charge screening in the double layer.(A) Provides the ion charge screening factor, f, as a function of distance, x, from the graphite surface.(B) Provides the electric potential, j, calculated from ion charge distributions in the same region.The colour scale from blue to red in (A) and (B) indicates increasing bulk ion concentrations across the entire concentration range: 0.2-9.2M, as listed in Table S2 of ESI.† (C) Provides the potential difference across the interface region (i.e. the PZC) as a function of c b NaCl .The simulation data, taken from the difference in j i in panel (B), are shown by the blue circles (the left y axis scale apply to these data only), and measurements from electrochemical experiments are provided by the black circles (the right y axis scale applies here).Statistical uncertainties of one standard deviation in the data are shown by the error bars.The dashed lines provide a moving average of the data and the solid lines are a linear fit to the data when c b NaCl < 3 M, the gradients for which are provided.
decreases (due to increased levels of screening in a less diffuse counter-ion charge cloud) and, combined with increasing c b NaCl , the relatively large negative gradient, dj PZC i / dc b NaCl , at the lowest concentrations reduces on increasing c b NaCl .At the highest concentrations, however, the ordering of ions leads to a small positive change in Ð (1 À f(x)) on increasing c b NaCl .This is sufficient to change the sign of dj PZC i /dc b NaCl .It follows that, for small increases to c b

Fig. 4
Fig. 4 Contributions to the electrochemical potential of electrolytes in the double layer region.(A) NaCl molalities (b b NaCl ) when the bulk mean molalities were as shown in the legend in units of m (mol kg À1 ).(B) The contribution to the mean ion molal electrochemical potential (m) from the first term in eqn (6).(C) The contribution of the second term in eqn (6) to m.The curve for the case when b b NaCl ¼ 3.5 mol kg À1 is highlighted by the bold line, and the shaded region indicates the position of the maximum peak in b b NaCl .

Fig. 5
Fig. 5 Ion assembly and diffusion in the double layer.(A) Average Na-Cl first-sphere coordination number (N Na-Cl ) as a function of x calculated from a CmMD simulation where c b NaCl ¼ 9.2 M. The black and red dashed lines indicate the maximum first Na and first two Cl densities in the concentration profiles for Na and Cl highlighted in Fig. 1.Inset is a configuration of a particularly large ionic network identified within the red dashed lines in x.The image is projected onto the yz plane; blue and cyan spheres represent Na + and Cl À ions and the green lines highlight ions that are directly coordinated.(B) Ion clusters observed in the bulk and double layer regions of simulations.The blue / red colour scale in the main plot indicates increasing c b NaCl from 1.2 to 9.2 M. The surrounding panels show purple and green probability densities in the same 2D space for clusters in the bulk and in the double layer, respectively, at the highest bulk concentrations shown by the legend.(C) Diffusion coefficients, D, for ions and water in 0.4 nm regions at the interface and in the bulk region of CmMD simulations (indicated by the subscripts i and b in the legends).The top panel provides the x component of D. Data have been scaled by 1 Â 10 À5 cm 2 s À1 and error bars show uncertainties of one standard deviation.

Fig. 6
Fig. 6 The effects of applied surface charges on the double layer when c b NaCl $ 1 M. (A and B) Provide ion concentrations as a function of distance from the electrode (positioned at x ¼ 0) measured in CmMD simulations when jsj ¼ AE0.92e nm À2 , where s is the surface charge density.(C) Provides the difference in the maximum ion concentrations (Dc) at the electrode as a function of s.Equivalent concentration differences, from distributions when no surface charge is applied (Dc ref ), were first subtracted, and shaded areas show uncertainties in the MD data due to changes using a smoothing window of 0.2 AE 0.1 nm.(D) Shows the integrals of 1 À f 0 , where f 0 is the electrode screening factor, calculated for all atoms in the double layer region in solutions.(E) Provides the capacitance (C) as a function of the potential difference across the double layer after subtracting the potential of zero charge (Dj 0 ).Experimental data are calculated from the integral of the differential capacitance (C d ) data in panel F for the 1 M case.(F) Provides the measured C d from experiments and the evaluated double layer capacitance from simulations (C d dl ) as a function of Dj 0 .The dashed blue line is a moving average over three data points (in blue) at positive/negative Dj 0 .This is shown more clearly in the inset of (F) over the range of potential differences used to measure C d (also inset and shown by the black data points).Where shown, error bars highlight the standard error of the mean from multiple trajectory window analyses (otherwise statistical uncertainties are of the size or smaller than the data points), with the grey region in (F) showing the uncertainty in the mean C d from three repeat experiments.In (E) and (F), data have been truncated for ease of comparison between simulations and experiments.