Alexandre Cauquoin
email: alexandre.cauquoin@awi.de

Camille Risi

Etienne Vignon

Étienne Vignon

Importance of the advection scheme for the simulation of water 1 isotopes over Antarctica by atmospheric general circulation models: a case study for present-day and Last Glacial Maximum with LMDZ-iso

Keywords: water stable isotopes, Antarctica, AGCM, advection, isotope-temperature gradient

Atmospheric general circulation models (AGCMs) are known to have a warm and isotopically enriched bias over Antarctica. We test here the hypothesis that these biases are consequences of a too diffusive advection.

Exploiting the LMDZ-iso model, we show that a less diffusive representation of the advection, especially on the horizontal, is very important to reduce the bias in the isotopic contents of precipitation above this area. The choice of an appropriate representation of the advection is thus essential when using GCMs for paleoclimate applications based on polar water isotopes. Too much diffusive mixing along the poleward transport leads to overestimated isotopic contents in water vapor because dehydration by mixing follows a more enriched path than dehydration by Rayleigh distillation. The near-air surface temperature is also influenced, in a more minor way, by the diffusive properties of the advection scheme directly via the advection of the air and indirectly via the radiative effects of changes in high cloud fraction and water vapor. A too diffusive horizontal advection increases the temperature and so also contributes to enrich the isotopic contents of water vapor over Antarctica through a reduction of the distillation. The temporal relationship, from Last Glacial Maximum (LGM) to presentday conditions, between the mean annual near-air surface temperature and the water isotopic contents of precipitation for a specific location can also be impacted, with significant consequences on the paleotemperature reconstruction from observed changes in water isotopes.

Introduction

Water stable isotopologues (hereafter designated by the term "water isotopes"), are integrated tracers of the water cycle. Especially, the isotopic composition recorded in polar ice cores enabled the reconstruction of past temperature variations (Jouzel, 2013 and references therein). For example, low accumulation sites that are typical on the East Antarctic Plateau (< 10 cm water-equivalent yr -1 ) provided the longest ice core records, allowing to reconstruct past climate over several glacial-interglacial cycles [START_REF] Jouzel | Orbital and millenial antarctic climate variability over the last 800,000 years[END_REF]. However, the interpretation of isotope signals remains challenging because of the numerous and complex processes involved (water vapor transport, fractionation during the phase changes in the water cycle, distillation effect…). This is particularly the case for Antarctica because this part of the world is subject to extreme weather conditions.

To improve our knowledge on the mechanisms controlling the water isotopes distribution, atmospheric general circulation models (AGCMs) enhanced by the capability to explicitly simulate the hydrological cycle of the water isotopes (H2 16 O, HDO, H2 17 O, H2 18 O) are now frequently used [START_REF] Joussaume | A general circulation model of water isotope cycles in the atmosphere[END_REF]Risi et al., 2010a;[START_REF] Werner | Stable water isotopes in the ECHAM5 general circulation model: Toward high-resolution isotope modeling on a global scale[END_REF]. Water isotopes in climate models have been used, for example, to better understand how the climatic signal is recorded by isotopes in polar ice cores at paleoclimatic time scales [START_REF] Werner | Isotopic composition and origin of polar precipitation in present and glacial climate simulations[END_REF].

However, some issues remain concerning the simulation of the climate over the Antarctic continent by AGCMs.

For example, they frequently present a near-surface warm bias over this area [START_REF] Masson-Delmotte | Past and future polar amplification of climate change: climate model intercomparisons and ice-core constraints[END_REF] and isotopic values in precipitation that are not depleted enough compared to observations [START_REF] Lee | Analysis of the global distribution of water isotopes using the NCAR atmospheric general circulation model[END_REF]Risi et al., 2010a;[START_REF] Werner | Stable water isotopes in the ECHAM5 general circulation model: Toward high-resolution isotope modeling on a global scale[END_REF]. This raises the question why many of the AGCMs have these warm and enriched in heavy water isotopes biases over Antarctica.

In this paper, we hypothesize that one part of these biases is associated with an excessively diffusive water vapor transport, i.e. transport that is associated with too much mixing. According to previous studies, the diffusive properties of the advection scheme in the AGCMs, on the horizontal as well on the vertical, can have an impact on the simulation of humidity and of its water isotope contents. On the horizontal, dehydration of air masses by mixing with a drier air mass leads to more enriched water vapor than dehydration by condensation and associated Rayleigh distillation [START_REF] Galewsky | An advection-condensation model for subtropical water vapor isotopic ratios[END_REF]. For the same reason, poleward water vapor transport by eddies (which act as mixing) leads to more enriched water vapor in Antarctica than transport by steady advection [START_REF] Hendricks | Space and time variation of d 18 O and dD in precipitation: Can paleotemperature be estimated from ice cores?[END_REF]. On the vertical, the excessive diffusion during water vapor transport seems to be the cause of the moist bias found in most AGCMs in the tropical and subtropical mid and upper troposphere, and of the poor simulation of isotopic seasonality in the subtropics [START_REF] Risi | Process-evaluation of tropospheric humidity simulated by general circulation models using water vapor isotopic observations: 2. Using isotopic diagnostics to understand the mid and upper tropospheric moist bias in the tropics and subtropics[END_REF]. The diffusivity of the advection scheme in the vertical has also important consequences on modeling of tracers like tritium by affecting greatly its residence time in the stratosphere, and so its downward transport from the stratosphere to the troposphere [START_REF] Cauquoin | Modeling the global bomb tritium transient signal with the AGCM LMDZ-iso: A method to evaluate aspects of the hydrological cycle[END_REF].

The goal of this paper is to test whether the warm and enriched biases in Antarctica are associated with an excessively diffusive water vapor transport, both on the horizontal and on the vertical. The diffusive character of the advection can be varied by modifying either the advection scheme or the resolution of the simulation, and we test both possibilities. Finally, we explore if a too diffusive water vapor transport can affect the temporal water isotopes -temperature slope between the Last Glacial Maximum (LGM, 21 ka) and present-day periods.

Model, simulations and data

Model and simulations

We use here the isotopic AGCM LMDZ-iso (Risi et al., 2010a) at a standard latitude-longitude R96 grid resolution (2.5° × 3.75°), and with 39 layers in the vertical spread in a way to ensure a realistic description of the stratosphere and of the Brewer-Dobson circulation [START_REF] Lott | The stratospheric version of LMDz: Dynamical climatologies, Arctic oscillation, and impact on the surface climate[END_REF]. Water isotopes are implemented in a way similar to other state-of-the-art isotope-enabled AGCM (Risi et al., 2010a). The isotopic composition of glacier Rglacier is calculated prognostically in the model. It is a precipitation-weighted average of the previous snow fall.

At each time step and in each grid box, it is updated as:

𝑅 glacier (𝑡 + 𝑑𝑡) = / glacier × 2 glacier (3) 4 567 snowfall × =3 / >?@ABCD 4 E F G snowfall × =3 , (1)
with hglacier = 20 kg.m -2 the height scale for the glacier, isosnowfall and H2Osnowfall the snowfall of isotopes and standard water in kg.m -2 .s -1 and dt = 30´60 s. No fractionation is assumed during the runoff and sublimation of glaciers ice. The model has been validated at global scale for the simulation of both atmospheric [START_REF] Hourdin | The LMDZ4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection[END_REF] and isotopic (Risi et al., 2010a) variables, and has been extensively compared to various isotopic measurements in polar regions [START_REF] Casado | Impact of precipitation intermittency on NAO-temperature signals in proxy records[END_REF][START_REF] Steen-Larsen | Continuous monitoring of summer surface water vapor isotopic composition above the Greenland Ice Sheet[END_REF], 2014, 2017;[START_REF] Bonne | The isotopic composition of water vapour and precipitation in Ivittuut, southern Greenland[END_REF][START_REF] Bonne | The summer 2012 Greenland heat wave: In situ and remote sensing observations of water vapor isotopic composition during an atmospheric river event[END_REF][START_REF] Touzeau | Acquisition of isotopic composition for surface snow in East Antarctica and the links to climatic parameters[END_REF][START_REF] Ritter | Isotopic exchange on the diurnal scale between near-surface snow and lower atmospheric water vapor at Kohnen station, East Antarctica[END_REF][START_REF] Stenni | Three-year monitoring of stable isotopes of precipitation at Concordia Station, East Antarctica[END_REF]. LMDZ-iso is also able to simulate the H2 17 O distribution [START_REF] Risi | Can we determine what controls the spatio-temporal distribution of d-excess and 17 O-excess in precipitation using the LMDZ general circulation model?[END_REF] but we do not consider it here because the limitations inherent to the AGCMs lead to strong uncertainties and numerical errors on the spatio-temporal distribution of this isotope. 2 years of spinup have been performed for all the simulations presented hereafter.

To quantify the effects of the prescribed advection scheme on water stable isotope values over Antarctica, we first performed three sensitivity simulations with LMDZ-iso under present-day conditions for the post spin-up period 1990-2008 (i.e. 19 model years), following the model setup from [START_REF] Cauquoin | Modeling the global bomb tritium transient signal with the AGCM LMDZ-iso: A method to evaluate aspects of the hydrological cycle[END_REF] i.e. simulations follow the AMIP protocol [START_REF] Gates | AMIP: The Atmospheric Model Intercomparison Project[END_REF], forced by monthly observed sea-surface temperatures and nudged by the horizontal winds from 20CR reanalyses [START_REF] Compo | The Twentieth Century Reanalysis Project[END_REF]: (1) one control simulation with the van Leer (1977) advection scheme (called VL), which is a second order monotonic finite volume scheme prescribed by default in the standard version of the model (Risi et al., 2010a); and two other simulations whose the van Leer advection scheme has been replaced by a single upstream scheme [START_REF] Godunov | Finite-difference methods for the numerical computations of equations of gas dynamics[END_REF] on (2) the horizontal plane (UP_xy) and on (3) the vertical direction (UP_z). Depending on one tunable parameter, the LMDZ model can be used with these 2 versions of the advection scheme according to the object of study [START_REF] Risi | Process-evaluation of tropospheric humidity simulated by general circulation models using water vapor isotopic observations: 2. Using isotopic diagnostics to understand the mid and upper tropospheric moist bias in the tropics and subtropics[END_REF]. The advection scheme in the simulations presented in the LMDZ-iso reference paper from Risi et al. (2010a) was set erroneously to the simple upstream scheme rather than to the van Leer's scheme (Risi et al., 2010b), and has little influence on their simulated spatial and temporal distributions of water isotopes at a global scale. However, as we will show here, this has considerable effect on the spatial distribution of these proxies over region with extreme weather conditions such as Antarctica. The 2-year spin-up time is enough to reach equilibrium. In the VL simulation, the globally and annually average values of temperature and d 18 O in precipitation for 1990 are 13.06°C and -7.53‰ respectively, very close to the average values over the whole period 1990-2008 (13.14°C and -7.57 ‰) within the average interannual variability of 0.15°C and 0.06 ‰. The conclusion is the same if we focus on the 60°S-90°S area instead: the average values of temperature and d 18 O in precipitation for the year 1990 are -17.89°C and -25.45 ‰ respectively, very close to the average values over the whole period 1990-2008 (-17.71°C and -25.34 ‰) within the average interannual variability of 0.43°C and 0.29 ‰.

The detailed description of the mixing ratio by the van [START_REF] Van Leer | Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection[END_REF] scheme and its comparison with the upstream scheme [START_REF] Godunov | Finite-difference methods for the numerical computations of equations of gas dynamics[END_REF] can be found in the Appendix A of [START_REF] Cauquoin | Modeling the global bomb tritium transient signal with the AGCM LMDZ-iso: A method to evaluate aspects of the hydrological cycle[END_REF]. To resume, the mixing ratio at the left boundary of box i, qi-1/2, is calculated as a linear combination of the mixing ratio in the boxes i-1 and i in the van Leer's scheme whereas in the upstream scheme qi-1/2 = qi-1. This means that in the upstream scheme, even if the air mass flux from grid box i-1 to grid box i is very small, the air that is advected into box i has the same water vapor mixing ratio as grid box i-1 as a whole. This makes the upstream scheme much more diffusive. An example of the effect of these 2 different advection schemes in a very idealized case is given in the figure 3 (top panel) of [START_REF] Hourdin | The Use of Finite-Volume Methods for Atmospheric Advection of Trace Species. Part I: Test of Various Formulations in a General Circulation Model[END_REF]. In a single dimension case, a tracer distribution is initially rectangle and is advected with a constant velocity such that Courant number is equal to 0.2. After 70 iterations, the initial rectangle shape is almost unchanged with the van Leer scheme, whereas it spreads into a flat Gaussian shape with the upstream scheme. Quantitatively, almost half of the tracer mass becomes outside the initial rectangle shape with the upstream scheme.

Increasing the grid resolution is equivalent to using an advection scheme that is less diffusive. Indeed, these finite-difference schemes are discretization methods and so depend on the chosen spatial resolution. To check that our findings and conclusions are consistent, we performed two more UP_xy and VL simulations under present-day conditions with the same configuration as presented above but at the R144 resolution (latitudelongitude grid resolution of 1.27° × 2.5°).

Finally, we evaluate the impacts of applying advection schemes with different diffusive properties on the temporal relationship between water isotopic contents of precipitation and mean air surface temperature, essential for paleo-temperature reconstructions. For this, we effected 12 years long (the two first years being used for the spin-up) UP_xy and VL simulations for present-day "free" (i.e. not nudged by the 20CR reanalyses, still following the AMIP protocol) and LGM conditions at the R96 resolution. For the LGM simulations, the PMIP3 protocol is applied [START_REF] Braconnot | Evaluation of climate models using palaeoclimatic data[END_REF]. Orbital parameters and greenhouse gas concentrations are set to their LGM values. ICE-5G ice sheet conditions are applied [START_REF] Peltier | Ice age paleotopography[END_REF]. LMDZ-iso is forced by the climatological sea-surface temperatures (SSTs) and sea ice from the IPSL-CM4 model [START_REF] Marti | The new IPSL climate system model: IPSL-CM4[END_REF]. The SST simulated by IPSL-CM4 for pre-industrial conditions (PI) has a global bias of -0.95 Kelvin with a cold bias in the mid-latitudes, a warm bias on the eastern side of the tropical oceans and on the Southern Ocean, and a particularly strong cold bias in the North Atlantic [START_REF] Hourdin | Impact of the LMDZ atmospheric grid configuration on the climate and sensitivity of the IPSL-CM5A coupled model[END_REF]. To avoid confusing these biases with LGM -present-day signals,

Observations

For analyzing the model performance over Antarctica under present-day conditions, we make use of the 1a) due to the very low temperatures. Because of the extreme cold and dry conditions at this area, one can see that the main disagreements between model outputs and observations are located at this place (Figure 1 and blue background of Table 1). Again, the isotopic outputs from the VL simulation are in better agreement with the observations (Table 1, blue background). These first results confirm that an excessively diffusive water vapor transport influences significantly the simulated isotopic and temperature values over Antarctica.

Temperature

The bias in temperature is deteriorated about in the same way when applying a more diffusive advection on the vertical direction or on the horizontal plane, as shown with the RMSE of annual mean temperature of 7.50, 7.31 and 6.60°C for the UP_z, UP_xy and VL simulations respectively (Table 1, red background). This tendency is the same when focusing on the East-Antarctic plateau. However, in average over the East Antarctic plateau, the temperatures of -30.51, -30.69 and -31.54°C from the UP_z, UP_xy and VL simulations are all within the spatial average of standard deviations of long-time temperature values (0.87°C). These values are all much warmer than the average observed temperature (-36.93°C). This shows that other factors than advection are responsible for the warm biais.

It has been suggested that the Antarctic warm bias in AGCMs could be linked to the general poor representation of the polar atmospheric boundary layer and related atmospheric inversion temperatures in these models [START_REF] Krinner | Studies of the Antarctic climate with a stretched-grid general circulation model[END_REF]. [START_REF] Cesana | How well do climate models simulate cloud vertical structure? A comparison between CALIPSO-GOCCP satellite observations and CMIP5 models[END_REF] have also shown that CMIP5 models generate too many highlevel clouds (i.e. above an altitude of 6.72 km), that can partly explain the overestimation of temperatures in Antarctica, due to their effect on downwelling longwave radiation. To go further, we compare the seasonal signals of surface temperature, near-surface thermal inversion (defined here as the difference between the 10mtemperature and the surface temperature), downward longwave radiative flux at the surface (LWdn), integrated water vapor column and high cloud fraction at EPICA Dome C (EDC: [75.10° S; 123.35° E]) from VL and UP_xy simulations with meteorological observations (see Figure 2 and the description of the data in section 2.2).

Both VL and UP_xy simulations have a warm bias at the surface (Figure 2a) and underestimate the near-surface inversion (Figure 2b), especially during the austral winter. This is mostly explained by an overly active turbulent mixing within stable boundary layers in LMDZ5 [START_REF] Vignon | Antarctic boundary layer parametrization in a general circulation model: 1-D simulations facing summer observations at Dome C[END_REF]. The disagreement with the observed nearsurface inversion at EDC is exacerbated in the UP_xy simulation because the modeled surface temperature is more affected by the change of advection scheme than the modeled 10m-temperature. This is consistent with a higher total LWdn flux over EDC in the UP_xy simulation than in the VL one (Figure 2c), by around 5 W.m -2 in summer and 8 W.m -2 in austral winter. The difference of LWdn between our two simulations can be due to two aspects: the fraction of high cloud and the amount of water vapor (e.g., [START_REF] Vignon | Modeling the Dynamics of the Atmospheric Boundary Layer Over the Antarctic Plateau With a General Circulation Model[END_REF]. The modeled highlevel cloud fraction over EDC for the period 2007-2008 is overestimated compared to the CALIPSO-GOCCP observations (Figure 2f). This finding is consistent with the results from [START_REF] Cesana | How well do climate models simulate cloud vertical structure? A comparison between CALIPSO-GOCCP satellite observations and CMIP5 models[END_REF] for CMIP5 models and with [START_REF] Lacour | How well are clouds simulated over Greenland in climate models? Consequences for the surface cloud radiative effect over the ice sheet[END_REF] for dry areas over Greenland. This result also concurs with an overestimated relative humidity with respect to ice in the troposphere compared to radiosoundings (not shown) during winter seasons, especially in the UP_xy simulation. The disagreement of high cloud cover is clearly enhanced in the UP_xy case, which is consistent with stronger LWdn in the UP_xy simulation. The comparison of LWdn under clear sky in our two simulations (Figure 2d) can give us information about the contribution of cloud and water vapor on the variations of total LWdn (Figure 2c). The stronger simulated LWdn in the UP_xy simulation is due mainly to the high cloud in winter (70 %) and to both high cloud and water vapor during summer (42.5 and 57.5 %, respectively). The effect of water vapor on downward longwave radiative flux is also confirmed by the higher amount of precipitable water over EDC in the UP_xy simulation (Figure 2e). From these results, we deduce that the warmer surface temperature in the UP_xy simulation is due to the higher air temperature (direct effect of the horizontal upstream advection) and to the radiative amplification from high clouds and in a lesser extent from water vapor (indirect effect of the horizontal upstream advection).

The near-surface warm bias in LMDZ-iso, which is most pronounced for the coldest temperatures (see Figure 3), has the consequence that the distillation is not strong enough. Some microphysical processes and kinetic fractionation at very low temperature can be missed too. These different aspects could contribute to an overestimation of the d 18 O and dD in precipitation over Antarctica. Finally, the 20CR reanalysis assimilate only surface observations of air pressure and use the observed monthly sea surface temperature and sea ice concentration as lower boundary conditions. These less strong constraints, compared to other reanalyses, may cause biases on the surface temperature over the poles (A. Orsi, personal communication), that can impact our isotopic delta values.

Spatial d 18 O-temperature relationship

We compare now our simulated spatial d 18 O-temperature relationship and d 18 O values for a given temperature to the ones from the data compiled by [START_REF] Masson-Delmotte | A review of Antarctic surface snow isotopic composition: observations, atmospheric circulation and isotopic modelling[END_REF]. Over the full temperature range, the spatial gradient is 0.83 ‰.°C -1 in the VL simulation, very close of the observed one (0.80 ‰.°C -1 ). We make the same comparison but by restricting the dataset to the ones on the East-Antarctic plateau (Figure 3). As noticed previously, the average modeled temperature over Antarctica is overestimated whatever the simulation considered. Especially, no simulated temperature reaches a value below -50°C. Yet, simulated d 18 O values are too depleted for the temperature range between -49°C and -43.5°C. As a consequence, a steeper modeled d 18 Otemperature gradient is observed for the lowest temperatures, giving a modeled global gradient of 1.24 ‰.°C -1 (thin orange line). If we restrict the fit to simulated temperatures higher than -43.5°C (vertical dashed line in

Figure 3), corresponding to the change of slope in the simulated d 18 O-temperature relationship (thick orange line), the simulated gradient (0.96 ‰.°C -1 ) is in more reasonable agreement with the one from the observations (0.85 ‰.°C -1 ).

We now discuss possible reasons that could explain why the simulated d 18 O-temperature slope is too steep at very low temperatures. First, it could be related to missing representation of fractionation during sublimation from the surface. As for all AGCMs equipped with water isotopes, fractionation at sublimation is not taken into account in LMDZ-iso. However, this effect would lead to further decrease of the water vapor d 18 O in polar region and hence contribute to an even steeper d 18 O-temperature slope at low temperature (hence further accentuate the mismatch). Second, the slope mismatch could be related to poorly represented kinetic fractionation. As in all the other models equipped with water isotopes, the parameterization of kinetic effect during vapor-to-solid condensation is represented empirically using a linear relationship between the supersaturation and the condensation temperature (Risi et al., 2010a). A modification of the temperature can thus induce some change in the d 18 O of the condensate, but this effect is of second order compared to the distillation effect explaining much of the slope between d 18 O and surface temperature. Third, the slope mismatch could be related to a poor representation of the atmospheric boundary layer and of its related inversion temperature [START_REF] Krinner | Studies of the Antarctic climate with a stretched-grid general circulation model[END_REF][START_REF] Masson-Delmotte | Past and future polar amplification of climate change: climate model intercomparisons and ice-core constraints[END_REF], as shown in the Figure 2. In LMDZ-iso, the warm bias in the simulated condensation temperature is smaller than that in the simulated surface temperature. Therefore, the water vapor masses continue to be distillated when moving away from the coast, while the cooling simulated at the surface from the coast to the remote region of the East Antarctic plateau is much less steep than in the reality.

Comparison of the different simulations

Effects of the diffusive properties of the advection scheme

We compare here the results from our different present-day simulations at a R96 grid resolution. The UP_z simulation (upstream vertical advection, Figure 1b) increases the bias a little in d 18 O, but its results stay relatively close of the d 18 O values from the VL simulation, indicated by the similar average values that differ only by 0.89‰ for all Antarctica (Table 1) that is smaller than the mean of the 60° S -90° S standard deviations. On the other hand, the d 18 O outputs from the UP_xy simulation (upstream horizontal advection, Figure 1c) display greater differences with the VL simulation ones, and so with the isotopic data, as revealed by the mean UP_xy -VL difference in d 18 O of 4.31‰. This is even more significant when focusing on the East-Antarctic plateau, with a respectively. The difference between the VL and UP_xy gradient shows an effect of diffusive properties of the large-scale transport on the distillation process. This difference between the modeled d 18 O-temperature gradients is reduced if we restrict to the temperature range above -43.5°C over the East-Antarctic plateau, with gradients of 0.93 and 0.96 ‰.°C -1 according to the UP_xy and VL simulations respectively (Figure 3, red and orange thick lines).

Since the simulated temperature difference between the UP_xy and VL simulations is in the error margin, i.e.

less than 1°C, we do not expect that the temperature difference explains the difference in spatial slopes. Rather, the larger slope in simulation VL is due to the larger relative contribution of Rayleigh distillation compared to mixing.

Effects of the horizontal grid resolution

We test now the hypothesis that to increase the horizontal resolution is equivalent to using an advection scheme that is less diffusive. The Antarctica-mean results are summarized in the Table 2. Compared to the UP_xy R96 simulation, the average value of d 18 O in precipitation is decreased by 4.31‰ when the advection scheme is improved (VL R96), and by 2.42‰ when the horizontal resolution is increased (UP_xy R144), in better agreement with the observations. The picture is the same for the dD outputs. The decrease of the mean modeled temperature values, smaller than the mean of the long-time standard deviations on the 60° S -90° S area, is the same by changing the advection scheme or by increasing the resolution: by 0.85°C and 0.76°C respectively. The best results are reached by improving both the advection scheme and the horizontal resolution at the same time, with model-data differences in temperature and d 18 O of 4.83°C and 0.88‰ respectively. This confirms that an increase of the horizontal resolution plays the same role as an improvement in the representation of the advection scheme on its horizontal plane. It is worth mentioning that the improvement in model-data agreement using a higher horizontal grid resolution is probably not only due to an improved representation of the advection in the model but, among others, to a better resolved Antarctic topography and near surface circulation. Our results are consistent with the study of [START_REF] Werner | Stable water isotopes in the ECHAM5 general circulation model: Toward high-resolution isotope modeling on a global scale[END_REF] that shows, for an increased horizontal resolution, a better agreement of the simulated isotopic delta values and of the water isotope-temperature gradient with the observations. The seasonality in modeled Antarctic temperature and precipitation, which influences the precipitation-weighted isotope values, is not altered by the change of advection scheme or horizontal grid resolution.

Effects on the LGM to present-day change in temperature

We evaluate now the effects of the diffusive properties of the advection scheme on the d 18 O -temperature temporal slope at different Antarctic locations between LGM and present-day. For this, we compare the VL and

UP_xy LGM simulations with their present-day "free" counterparts (i.e. without nudging). Figure 4a shows VL simulated temporal slopes in d 18 O -temperature at each location as Dd 18 O/DT. Over East Antarctica, the simulated temporal slope becomes larger inland, going from near to 0 ‰.°C -1 on the coast to 2 ‰.°C -1 on the deep East-Antarctic plateau. Over West Antarctica, the temporal slopes are more heterogeneous and show variations from 0.1 to 1 ‰.°C -1 . The slopes at WDC (Wais Divide Core, square symbol), Siple Dome (triangle symbol), Vostok (circle symbol) and EDC (cross symbol) are of 0.80, 0.66, 0.73 and 0.37 ‰.°C -1 respectively (calculated by considering the 9 grid cells centered on each drill location). Figure 4b shows the differences between the UP_xy and VL simulated slopes. As the spatial slopes from the VL simulation are higher than the UP_xy one, it is expected that the LGM to present-day temporal slopes are larger in the VL simulation than in the UP_xy one because of the colder temperatures and of the higher relative contribution of the Rayleigh distillation. This is the case on the western part of the continent where the UP_xy temporal slopes are smaller than the VL ones by up to 0.6‰.°C -1 .

Conclusions

We have tested with LMDZ-iso if the warm and isotopically enriched biases in Antarctica, frequently observed in the AGCMs, are associated with the diffusive property of the advection scheme. The simulated water isotope contents in Antarctica are very sensitive to the diffusive character of the water vapor transport on the horizontal plane. The higher the contribution of mixing (i.e. diffusion), the more enriched the precipitation. These findings are even more striking for the East-Antarctic plateau where the main ice cores allowing paleoclimate reconstructions are located. Moreover, because the diffusive character of the large-scale transport influences the temperature in this region, even in a light way, this has an impact on the modeled water isotopic composition through the Rayleigh distillation. So, we conclude here that the excessive numerical diffusion has a large influence on the enriched isotopic bias. For the spatial isotope-temperature relationship over the East-Antarctic plateau observed in LMDZ-iso, this latter is improved for the temperatures above -43.5°C, in more reasonable agreement with the observations. At the lowest temperatures (i.e still over the East-Antarctic plateau), that the model is not able to reach, the non-linearity observed in our simulations (the spatial d 18 O -temperature relationship is steeper for the lowest temperatures) can be unlikely explained at first order to missing or poorly represented kinetic fractionation. One can speculate that the water masses continue to be distillated when moving away from the coast, hence depleting the water vapor in heavy isotopes while the modeled temperature decrease from the coast to the remote region of the East Antarctic plateau is much less steep than in the reality.

This more pronounced effect in the case of a more diffusive horizontal advection scheme could be due to the deteriorated representation of the inversion temperature. The temporal isotope-temperature relationship at some locations in Antarctica can be influenced by the diffusive properties of the advection scheme on its horizontal domain. As for the spatial gradient, an excessive numerical diffusion has the consequence to decrease the isotope-temperature temporal gradient, leading to a wrong estimation of the LGM to present-day temperature change deduced from observed d 18 O. Our study demonstrates that a representation of the advection scheme in the AGCMs taking into account water isotopes and isotopic gradients, especially on the horizontal domain, is an important step toward a more realistic modeling of water isotopes over Antarctica.

Another way to improve this aspect is to increase the spatial resolution, which has the same effect as applying a less diffusive advection scheme on the water isotopic composition and the temperature. This study shows again the importance of using water stable isotopes in GCMs for the evaluation and quantification of the processes influencing the hydrological cycle, including the advection of water vapor. We expect our main results (excessive diffusive advection leading to warmer temperatures, moister boundary layer and more enriched water vapor) to be robust and to hold in other models as well. However, the quantitative response to the advection scheme may be modulated by the representation of boundary layer processes and high-cloud microphysics in each model. LGM-PD

d 18 OUP_xy (‰)
LGM-PD

d

  we use the SSTs from an IPSL PI simulation in the following way to cancel out the biases in the IPSL model common to both the LGM and PI simulations (see details inRisi et al., 2010a): SST forcing = SSTLGM -SSTPI + SSTAMIP. We set the sea surface d 18 O to +1.2‰ and assume no glacial change in the mean deuterium excess in the ocean, as inRisi et al. (2010a). Note once again that the 2-year spin-up time is enough to reach equilibrium. In the VL present-day "free" simulation, the globally and annually averaged values of temperature and d 18 O in precipitation for the first post spin-up model year are 12.54°C and -8.27 ‰ respectively, very close to the average values over the 10 years of simulation(12.54°C and -8.34 ‰) in comparison to the interannual variabilities of 0.06°C and 0.06 ‰. The conclusion still holds if we focus on the 60°S-90°S area instead: the average values of temperature and d 18 O in precipitation for the first post spin-up year are -22.03°C and -26.08 ‰ respectively, very close to the average values over the 10 years of simulation (-22.04°C and -26.04 ‰) within the average interannual variability of 0.10°C and 0.21 ‰. We express the isotopic composition of difference water bodies in the usual d-notation as the deviation from the Vienna Standard Mean Ocean Water (V-SMOW). So for H2 18 O, the d 18 O value is calculated as d 18 O = ([H2 18 O]/[H2 16 O])sample / ([H2 18 O]/[H2 16 O])V-SMOW -1) × 1000. Long-time mean d values are then calculated as precipitation-weighted mean. For the quantitative model-data comparisons, we retrieve the model values at data geographical coordinates by bilinear interpolation. Without such an interpolation, i.e. considering the nearest grid point instead, the root-mean-squared errors (RMSE) of mean d 18 O and temperature from the VL simulation differ only by 0.5‰ and 0.14°C compared to the results presented below, so the uncertainty associated with the model-data co-location is small.

Figure 1

 1 schemes for atmospheric boundary layer, inversion temperatures, we use additional datasets for the EPICA Dome C station (EDC: [75.10° S; 123.35° E]): the surface temperature, the 10m-temperature and the downward longwave radiative flux at the surface (LWdn) over the period 2011-2018 thanks to the CALVA program and the Baseline Surface Radiation Network (BRSN)(Vignon et al., 2018 and references therein); and the observed precipitable water from radiosoundings data for the period 2010-2017. We are aware that the observations periods are not the same than our model period(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008), giving a possible bias in the model-data comparison. For the cloud cover, we use the CALIPSO-GOCCP observations[START_REF] Chepfer | The GCM Oriented CALIPSO Cloud Product (CALIPSO-GOCCP)[END_REF] over the period2007-2008. We compare these data with the high-level cloud fraction over EDC simulated by LMDZ-iso through the COSP package[START_REF] Bodas-Salcedo | COSP: Satellite simulation software for model assessment[END_REF], which allows to simulate what the CALIOP instrument on-board CALIPSO would measure if it was into orbit above the simulated atmosphere. The difference in LWdn between the simulations VL and UP_xy at EDC is also present all over the East Antarctic plateau, confirming that EDC is representative of this area.

  model-data difference in d 18 O reaching 20‰ at some locations. The annual mean d 18 O and dD values from the UP_xy simulation are increased by 6.39‰ and 44.26‰ compared to the VL simulation average values, in less agreement with the observations as shown by their respective RMSE values (Table1, blue background). It shows that the diffusive property of the advection scheme on the horizontal plane is essential to better model the water isotope distribution, especially over Antarctica. To go further, one can also compare the d 18 O values at a fixed temperature for the UP_xy and VL simulations (Figure3, red and orange crosses respectively). The d 18 O in precipitation for a temperature of -32°C over the East-Antarctic plateau is already smaller by 6.1‰ in the VL simulation. This very significant difference in initial d 18 O can be attributed to the proportion of mixing against distillation that affects the water vapor during its transport. This lends support to our hypothesis that too much diffusive mixing along the poleward transport leads to overestimated d 18 O because dehydration by mixing follows a more enriched path than dehydration by Rayleigh distillation[START_REF] Hendricks | Space and time variation of d 18 O and dD in precipitation: Can paleotemperature be estimated from ice cores?[END_REF][START_REF] Galewsky | An advection-condensation model for subtropical water vapor isotopic ratios[END_REF]. We expect the relative contribution of mixing vs. distillation to have the largest impact on d 18 O at latitudes where eddies are the most active. This is probably why the d 18 O difference between the VL and UP_xy simulation becomes large in mid-latitudes, over the austral ocean before arriving at the Antarctica coast (not shown), hence the difference in "initial" d 18 O in Figure3.As noticed in section 3.1.2, all our simulations overestimate the average temperature in Antarctica and even more on the East-Antarctic plateau. A more diffusive advection on the horizontal or on the vertical increases the mean temperature value by 0.85 and 1.03°C respectively compared to the VL result. To explain such an influence of the advection on the temperature over Antarctica, even secondary, one can hypothesize that the Antarctic continent is better isolated, and so colder, when the advection of the model is less diffusive. If we focus now on the link between the temperature and the d 18 O in precipitation over all the continent, the d 18 O-temperature gradients according to our different R96 simulations UP_z, UP_xy and VL are at 0.79, 0.69 and 0.83 ‰.°C -1

Figure 1 :

 1 Figure 1: Map of Antarctica showing (a) the observed δ 18 O values from the compilation by Masson-Delmotte et al. (2008), (b) the difference between the simulated δ 18 O in precipitation and the δ 18 O observations for the UP_z, (c) UP_xy and (d) VL simulations. The bold black line shows the contour of 2500 m above sea level elevation.

Figure 2 :

 2 Figure 2: Multi-year monthly mean variations of (a) surface temperature, (b) near-surface thermal inversion (defined as the difference between the 10m-temperature and the surface temperature), (c and d) total and clear-sky component of downward longwave radiative flux at the surface, (e) precipitable water and (f) high cloud fraction above EPICA Dome C. The red and blue curves correspond to the VL and UP_xy simulation results, respectively. A comparison with observations (see section 2.2) is made when possible (dashed black lines).

Figure 3 :

 3 Figure 3: Relationship between δ 18 O and temperature on the East-Antarctic plateau according to the observations (black) and the UP_z (blue), UP_xy (red) and VL (orange) simulations. For each simulation outputs, two linear regressions have been conducted: one on the full East-Antarctic plateau dataset (thin lines) and one on the same dataset without the temperatures below -43.5°C (bold lines) indicated by the vertical dashed line. The corresponding formulas of these latter are also shown.

Figure 4 :

 4 Figure 4: (a) Temporal slope (Dd 18 O/DT) between the present-day (PD) and LGM according to the VL simulations. (b) Difference with the temporal slope deduced from the UP_xy simulations. The grey areas indicate where the LGM-PD difference in d 18 O is lower than the standard deviation of the interannual variability. Square, triangle, circle and cross symbols represent the WDC, Siple Dome, Vostok and EDC sites respectively.

Table 3 :

 3 Sites of interest and their geographical coordinates, observed LGM to present-day changes in d 18 LGM to present-day changes in d 18 O (UP_xy and VL simulations), simulated temporal d 18 O-temperature slopes and VL-UP_xy difference in reconstructed temperature change deduced from the observed d 18 O changes and the simulated temporal slopes. This difference is calculated as -(slope_UP_xy -slope_VL) ´ Dd 18 Oobs / slope_VL 2 . Negative signs indicate that when using the UP_xy temporal slope, we overestimate the LGM cooling compared to using the VL temporal slope. The model values are based on the spatial averages over the 9 grid cells surrounding the ice cores geographical coordinates.

  ). This is confirmed by the smaller root-mean-squared error of modeled d 18 O in precipitation from the VL simulation, calculated as the difference between the observed annual mean values and the LMDZ-iso results (RMSE = 4.47‰, i.e. 12.2% of the observed mean Antarctic d 18 O value). The results from the VL simulation for the other isotopic variable dD is also the closest of the observations with a RMSE of 40.93‰ (Table1, red background). Our simulated d 18 O in precipitation is very sensitive to the choice of the advection scheme on the horizontal plane,

with more enriched values when a more diffusive advection scheme is applied (Figure

1c

). This is reflected by the mean d 18 O model value and RMSE, higher by 3.16 and 4.31 ‰ than the VL simulation values. On the contrary, the sensitivity of Antarctic d 18 O values to the diffusive properties of vertical advection is weak (Figure

1b

) with an RMSE very close of the VL one (4.84 ‰). The results from the VL simulation for the other isotopic variable dD is also the closest of the observations with a RMSE of 40.93‰ (Table

1

, red background). According to the observations, the East-Antarctic plateau is where the water isotope values are the lowest (mean d 18 O below -40‰, Figure

  On the East-Antarctic plateau, the VL d 18 O -temperature slopes are larger between the longitudes 100°E and 140°E than in UP_xy as well. The slope values are smaller in VL than in UP_xy only over the center of the continent, which corresponds to the lowest simulated temperature values and so where the gradient between the d 18 O and temperature values is much steeper (see 2 nd paragraph of section 3.1.3), especially for LGM conditions.We highlight now the consequences that an excessively diffusive horizontal advection may have on the LGM cooling reconstructed from simulated temporal slopes. For this, we calculated the difference in the LGM to present-day change in temperature deduced from the observed changes in d 18 O and simulated slopes (Dd 18 O/DT) by VL and UP_xy, at WDC and Siple Dome stations for West Antarctica, and Vostok and EDC stations for East Antarctica (Table3). For WDC, Siple Dome and Vostok stations, the average changes in d 18 O between LGM and present-day are larger in the VL simulation, -11.98, -9.88 and -3.90 ‰ respectively, than in the UP_xy simulation (-7.52, -6.24 and -2.07 ‰ respectively). For the EDC station, the difference in LGM to present-day change in d 18 O between our VL and UP_xy simulations can be considered in the margin error (-1.97 and -1.70 ‰ for VL and UP_xy simulations respectively). The values deduced from the VL simulation are also in better agreement with the observations (Table3). As noticed previously, the temporal d 18 O-temperature gradients at WDC, Siple Dome, Vostok and EDC are also larger in the VL simulation by 57%, 65%, 92% and 26% respectively, compared to the values deduced from the UP_xy simulation. As a consequence, if one applies the slope from the UP_xy simulation instead of the VL one, we would overestimate the present-day-to-LGM cooling at these stations by 4.34, 5.83, 3.76 and 2.95°C respectively. Even if the change in the temporal d 18 O-temperature gradient at EDC is relatively smaller than for the other stations, the consequence on the reconstructed present-day-to-LGM cooling value is still important. These results show that a more diffusive advection scheme on the horizontal plane can affect greatly the LGM to present-day temperature change deduced from observed d 18 O in precipitation.

Table 1 : Observed and simulated annual mean values of temperature (T), δ 18 O and δD for the full Antarctic dataset (red background) and restricted to the East-Antarctic plateau (blue background), and the corresponding RMSE.

 1

		Mean	Mean	RMSE	Mean	RMSE	Mean	RMSE
		data	UP_z	UP_z	UP_xy	UP_xy	VL	VL
	T (°C)	-36.93	-30.51	7.50	-30.69	7.31	-31.54	6.60
	d 18 O (‰)	-36.76	-34.85	4.84	-31.43	7.63	-35.74	4.47
	dD (‰)	-289.62 -272.28	43.76	-251.34	62.00	-279.49	40.93
	T (°C)	-47.46	-39.49	8.38	-39.88	7.99	-40.71	7.23
	d 18 O (‰)	-46.77	-42.27	5.03	-37.37	9.69	-43.76	3.80
	dD (‰)	-366.98 -325.37	43.79	-291.99	76.44	-336.25	33.58

Table 2 : Comparison of the observed annual mean values of temperature (T), δ 18 O and δD for the full Antarctic dataset with four different LMDZ-iso simulations, combining different horizontal resolutions (R96 and R144) and different advection schemes (UP_xy and VL).

 2

		Mean data	Mean UP_xy	Mean VL	Mean UP_xy	Mean VL
			R96	R96	R144	R144
	T (°C)	-36.93	-30.69	-31.54	-31.45	-32.10
	d 18 O (‰)	-36.76	-31.43	-35.74	-33.85	-37.64
	dD (‰)	-289.62	-251.34	-279.49	-267.17	-289.55

Acknowledgements. We thank A. Landais for her useful suggestions on this manuscript, J.-B. Madeleine and C. Genthon for their kind support about inversion temperatures over EDC, and R. Guzman, C. Listowski and H. Chepfer for their help with the GOCCP data. This work was granted access to the HPC resources of IDRIS under the allocation 0292 made by GENCI. The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/20072013)/ERC grant agreement no. 30604. We thank Christophe Genthon and the CALVA program for acquiring and distributing meteorological data at Dome C (http://www.lmd.jussieu.fr/~cgenthon/SiteCALVA/CalvaBackground.html) with the support of the french polar institute (IPEV). We also acknowledge the WCRP-BSRN network and Angelo Lupi for the dissemination of radiation data. Radiosoundings data at Dome C are freely distributed in the framework of the IPEV/PNRA Project "Routine Meteorological Observation at Station Concordia"-www.climantartide.it.


Thomas Dufils

Guillaume Jeanmairet

Benjamin Rotenberg

Michiel Sprik

Mathieu Salanne

Simulating electrochemical

Keywords:

Despite many advances over the past decades [1,2], the efficient simulation of full electrochemical cells at the molecular scale, using electronic structure based calculations, remains a daunting task. This is due to their slab structure, since the minimal experimental setup consists of an electrolyte between two electrodes. The system is generally simplified by simulating one interface only, but the main conceptual difficulty is to find a way to charge the electrode surface at fixed composition of the electrolyte. Several methods have recently emerged [3][4][5], where the system is allowed to exchange electrons with a reservoir at fixed voltage (grand-canonical approach), but they all rely on the use of continuum descriptions for the electrolyte. These models, which are generally based on a Poisson-Boltzmann theory [6], remain mostly qualitative and an atomistic description would be preferable (this is also true because the solvent may actively participate to electrochemical reactions [7,8]). This is almost impossible to do since it would be necessary to remove/insert ions to counterbalance the electrode charge fluctuations.

Here we propose an alternative route to simulate electrochemical cells. Our approach is based on the coupling of a finite field with a system consisting in an electrolyte and a single electrode. Finite fields methods, developed in the framework of the modern theory of polarization, consist in imposing a macroscopic field (electric field [9], polarization [10] or electric displacement [11]) via an extended Hamiltonian accounting for the interaction between the system and the fixed field. They have recently been adapted and applied to the study of electrical double layers at solid/liquid interfaces, and more precisely charged [12] or polar [13] insulators/electrolyte interfaces. Electrochemical systems are by nature more complex since they involve metallic electrodes whose charge distribution is not fixed but depends on the surrounding medium and the applied potential. Due to the long sim-ulation times related to the relaxation of the electrical double layer (which is typically longer than the nanosecond), we establish here a proof of concept by using a classical molecular dynamics (MD) setup. Indeed, even if metals can only be accurately described in the framework of quantum mechanics, models have been developed to reproduce the electrostatics in classical or mixed quantum/classical simulations [14,15]. In this contribution we extend the finite electric field method to an electrolyte interacting with such a model metallic electrode.

In the following, we focus on the description of the electrostatics, the Van der Waals interactions being represented by the conventional Lennard-Jones model. The charge density in any point of space is given by

ρ(r) = N i=1 q i δ(r -r i ) + M j=1 q j η 3 π 3/2 exp -η 2 (r -r j ) 2 (1)
where the first term is the contribution of the electrolyte, represented by a distribution of point charges with q i the partial charge of the atom i ∈ [1, N ] and r i its position; δ is the Dirac distribution. The second term represents the atoms of the metallic electrode, in which each site j ∈ [1, M ] is immobile (with position r j ) and carries a charge q j which is spatially distributed following a gaussian charge distribution of width η -1 . In order to represent the metallic character of the electrodes, the latter charges are allowed to fluctuate in response to the electrolyte fluctuations and thus are part of the microscopic degrees of freedom [14,16]. The Hamiltonian of the system is written as

H P BC = K(p N ) + U (r N , q M ), (2)
where K(p N ) is the kinetic energy which depends on the ion momenta p N = {p 1 ...p N } and U (r N , q M ) the potential energy, which depends on the ion positions r N = {r 1 ...r N } and charges of the electrode q M = {q 1 ...q M }.

We have appended a superscript PBC (periodic boundary conditions) to indicate that the electrostatic energies and forces are computed using standard Ewald summation to account. All the atoms within an electrode are held at constant potential by enforcing the following condition on each atom

Ψ P BC j = ∂U P BC C ∂q j = Ψ J (3)
where Ψ J is the prescribed potential of electrode J to which the atom j belongs and U P BC C is the Coulombic contribution to the energy, given by:

U P BC C = 1 2 ρ(r)ρ(r ) 4π 0 |r -r | drdr . ( 4
)
Formally, solving the set of self-consistent equations given by Equation 3 is equivalent to minimizing U P BC C j Ψ J q j with respect to the charges. Since this function is quadratic in the fluctuating charges q j (see Eq 4), the minimization can be efficiently performed with conjugate gradients. Note that we add an additional constraint by forcing the sum of the electrode charges to be null [17].

FIG. 1: Top: electrolyte-centered supercell (2D PBCs), with fixed applied potential between two electrodes (∆Ψ = Ψ right -Ψ left ). Bottom: conductor-centered supercell (3D PBCs), with a finite field E and a single electrode in which the potential of the atoms is set to the same value Ψ electrode .

The minimal components of an electrochemical cell are two electrodes and an electrolyte between them. In classical MD, the computational cost is not prohibitive so it is relatively easy to simulate complete systems instead of a single electrode surface. The conventional setup to simulate such systems is illustrated in the top panel of Figure 1. In the following we will refer to this system as the electrolyte-centered supercell (ECS). It is simulated using the 2D Ewald summation [18,19] since the two electrodes are held at different potentials (Ψ J = Ψ left or Ψ right ). We note x and y the two directions along which PBCs are used in this setup.

Finite field (E) simulations can be performed using the extended Hamiltonian introduced by Stengel and Vanderbilt in [11] and is written as:

H E = H P BC -Ω P • E, (5)
where H P BC is the Hamiltonian defined by Eq 2, Ω is the volume of the supercell and P the polarization per unit volume. In the modern theory of polarization, the dipole moment of a unit cell involving infinite periodic systems is viewed as a multivalued quantity, since it depends on the choice of the position of the periodic boundaries. Nevertheless, this is not a significant issue since only differences in polarization matter in the dynamics and in the calculation of physical properties, in practice via the itinerant polarization [20]:

P itinerant (t) = P itinerant (0) + 1 Ω N i=1 q i ∆r i (t), ( 6
)
where ∆r i is the displacement of the atom between time t = 0 and t for the "unfolded" trajectory, i.e. not taking jumps in position (hence polarization) across the periodic boundaries. The system must be periodic in the direction in which the finite electric field is applied, which implies the use of 3D PBCs. A field E corresponds to a drop of Poisson potential across the cell ∆Ψ = -EL z where L z is the length of the box in the direction of the field. From the practical point of view, a consequence of the use of 3D PBCs is that now the two electrodes of the ECS necessarily merge, yielding a single electrode at fixed potential Ψ J = Ψ electrode . The simulation cell can then be represented with the electrode at its center, yielding the conductor-centered supercell (CCS) shown in the bottom panel of Figure 1. Contrarily to the constant applied potential simulations, it is now the presence of the finite field that induces the polarization of the electrode and potential drop at the two electrode/electrolyte interfaces.

Coupling fluctuating charges to model conductors with finite field simulations requires two important adaptations of both methods. On the one hand, the cell polarization includes a contribution from the fluctuating charges as:

P conductor (t) = 1 Ω M j=1 q j (t)r j . (7)
This additional term does not depend on the position of the electrode inside the supercell since we enforce j q j = 0. On the other hand, the determination of the partial charges via the constraint of fixed potential includes an additional contribution to the electrostatic energy and potential due to the finite field. From the extended Hamiltonian (Eq 5) and the expression of the electrode contribution to polarization (Eq 7), one obtains the generalization of Eq 3 as:

Ψ j = Ψ P BC j -r j • E = Ψ J (8)
As an extension of the constant applied potential case, the self-consistent expressions given by Eq 8 are now solved by minimizing U P BC C -Ω P • Ej Ψ J q j with respect to the charges. Since P is a linear function of the charges, the performances of the conjugate gradient minimizer are not affected by this setup.

In order to test this new approach, we simulate two systems, one at constant applied potential between two distinct electrodes and the other with the finite field method and a single fixed-potential electrode as shown on Figure 1. The electrode(s) consist in a model structure made of a cubic crystal with the NaCl lattice constant, with the (111) plane facing the liquid. The intermolecular interactions consist in electrostatic interactions and Lennard-Jones potentials using Lorentz-Berthelot mixing. The electrode sites Lennard-Jones parameters are the ones of Cl - [21]. The electrolyte consists in an aqueous solution of NaCl composed of 603 SPC/E water molecules [22] and 20 ion pairs [21]. The cross-sectional area is 2.20 nm 2 and the length along the z axis L z = 11 nm.

The simulations were performed with a timestep of 2 fs in the NVT ensemble at 298 K using a Nosé-Hoover thermostat with a coupling constant of 0.4 ps. The systems were equilibrated during 10 ns before a production run of 10 ns.

For the finite field simulations, we use the CCS configuration with 3D PBCs. The single electrode, which is made of 12 planes of atoms, is set at null potential and the field E ranges from 0 to 2 V nm -1 . For the constant applied potential simulations, we setup the ECS configuration with 2D PBCs. In this setup the two electrodes are of equal dimensions (i.e. 6 planes of atoms each), kept under constant potentials Ψ left and Ψ right such that ∆Ψ = Ψ left -Ψ right = -EL z , using the same values for E as above. In both series of simulations, the value of the η parameter for the gaussian charges has been set to 0.5052 Å -1 following ref. [19]. Electroneutrality is enforced during the charge calculation process [17].

A first validation is provided by comparing the polarization of the electrodes with the two setups. To do this we split the single electrode in the CCS setup in two parts. This is easily made since one side accumulates positive charge and the other is exactly opposite, while the centre is almost neutral. Firstly, the average FIG. 3: Charge density within the electrode and the electrolyte for a finite field of 0.1 V nm -1 or the corresponding applied potential of 1.1 V. The contributions of water molecules and ions to the total charge distribution are also shown; they were shifted by -0.01 e Å-3 for a better readability. The grey-shaded areas correspond to the electrodes (the CCS was translated in order to match with the ECS setup).

accumulated charge Q on the positive side is compared with the ECS setup on Figure 2A for a wide range of applied fields (potentials). The agreement between the two methods is excellent since the points are almost superimposed. This is true not only for voltages up to 6 V where the charge increases linearly with the applied potential (which reflects a constant differential capacitance [16]), but also up to 20 V for which the variation is not linear anymore (see the inset). The probability distributions of the instantaneous values of this quantity (Figure 2B) are also identical, which shows that the two methods sample the same configuration space. Finally, Figure 2C illustrates the instantaneous electrode charges for a given electrolyte configuration. The two systems are indistinguishable, showing that even at the local scale the finite field method yields a correct representation of the electrode/electrolyte interface. When a finite field is applied to a bulk liquid, due to the PBCs the charges cannot accumulate in a specific region so that there is a net electric field in each point of the supercell. Here the presence of a blocking surface (the electrode) results in the formation of polarized layers on the electrolyte side. The polarization arises from two mechanisms: i) reorientation of the water molecules ii) local charge imbalance by accumulation of one ionic species and depletion of the other. The structure adopted by the liquid may be compared with the case of constant applied potential simulations. As shown on Figure 3, the agreement is again very good for the variation of the charge density and its splitting between water and ionic contributions across the simulation cell. This validates further the use of the single electrode in order to study electrochemical interfaces using 3D PBCs. Additional tests on the electric field and Poisson potential for several applied voltages are provided in the Supplementary Information; they all show the same level of accuracy.

In conclusion, we have demonstrated the possibility to simulate metal slabs with fluctuating surface charge and in the presence of an explicit electrolyte that counterbalances the charge. This is done through the combination of two methods: An applied finite field which polarizes the cell and a constant potential electrode that screens the field in bulk, leading to the formation of two independent electrochemical interfaces. The net charge on one electrode surface is the opposite of the net charge on the other, which maintains the overall charge neutrality of the metal slab. The electrode surface charge fluctuations are compensated by the adsorption of ions from the electrolyte. The method is validated through extensive comparisons with simulations using constant applied potential between two separated electrodes (and 2D PBCs). From the practical point of view, it is much easier to implement in classical MD packages since it avoids the introduction of 2D PBCs. From the computational point of view, the method is also more efficient: The simulation time is reduced by approximately 15 %. Last, but not least, it could much be more easily applied to the case of ab initio molecular dynamics, which would open the door for the first-principles simulation of electrochemical reactions occurring at electrodes, in the presence of an explicit electrolyte.

Supplementary Information

Constant potential and applied electric field

As seen in the main text, fixed field and constant potential agree very well on the charge distribution (Fig 4 A). From here one obtains the electric field profile in the z direction, averaged on the x and y directions, using the Maxwell-Gauss equation

dE z dz = 4πρ(z), (9)
with the condition E z = E, the applied field. Fig 4 B shows that both methods agree once again. In particular, they both display a null electric field in the bulk electrode and electrolyte as expected for conductors. This is highlighted when we calculate from the electric field profile the Poisson potential profile displayed on fig 4 C, since the potential is flat in the electrode for both methods. We also observe the same potential drop for both methods between the electrode and the electrolyte at each of the interfaces. Thus the capacitance of each half-electrode in the 3D case is the same as the 2D equivalent electrode.

The agreement between applied field and constant potential is thus extended to the Poisson potential.

Field profiles for the single elctrode

One feature we need to check is if the single electrode under an applied field does behave like a perfect metal. This implies that the electric field in the bulk of the electrode is really null for every applied electric field. The electric field profile is displayed for the single electrode for multiple values of the applied electric field on fig 5 .  Whatever the applied field, the field in the bulk electrode is equal to zero after the second atomic layer. A closer look on this field may be obtained by the calculation of the Poisson potential. This profile is displayed for the single electrode for multiple values of the applied electric field on fig 6. In the electrode, the Poisson potential profile is flat, corresponding indeed to a zero electric field in the elec- Electric field E z (V/nm)

E = 0.5 V/nm E = 0.4 V/nm E = 0.3 V/nm E = 0.2 V/nm E = 0.1 V/nm E = 0.0 V/nm
FIG. 5: Electric field profile for the single electrode under an applied electric field (corresponding to potential drop across the cell of 0 to 5.5 V). The gray-shaded area corresponds to the electrode. Poisson potential

Ψ (V) E = 0.5 V/nm E = 0.4 V/nm E = 0.3 V/nm E = 0.2 V/nm E = 0.1 V/nm E = 0.0 V/nm
FIG. 6: Poisson potential profile for the single electrode under an applied electric field (corresponding to potential drop across the cell of 0 to 5.5 V). The gray-shaded area corresponds to the electrode.

trode, which does behave like a perfect metal. The plateau observed allows to assign a value of the Poisson potential to the single electrode. We also notice that because of the charge density fluctuations in the electrolyte, the potential difference ∆Ψ between the two bulk part of the electrolyte on each side of the electrode is not strictly equal to -EL z . This slight difference is otherwise not significant and is lower than the fluctuations of ∆Ψ.

Fluctuation-dissipation relation

The capacitance of the system sudied in the main text may be obtained from the Q(∆Ψ) relation displayed in the main text, and gives a value of C = 1.94 ± 0.11 µF.cm -2 . It may also be calculated using the fluctuation-dissipation theorem:

C = 1 S Q 2 -Q 2 k B T , ( 10
)
where S is the cross-section and T the temperature. The obtained values and a comparison with the fit of the Q(∆Ψ) relation are displayed on figure 7, showing that the two approaches are equivalent. This is not surprising since we observed in the main text that the charge distributions for constant electric field and constant applied potential are very similar in terms of both average values and standard deviation. This confirms that the capacitance does not depend on the applied voltage in the range [0:5V], that we could infer from Q(∆Ψ) relation. Capacitance obtained using the fluctuationdissipation theorem for the finite electric field (red) and constant potential (blue). The grey-shaded area corresponds to the confidence interval obtained through the Q(∆Ψ) relation

1 BFIG. 2 :

 12 FIG.2:The single electrode behaves as two distinct constant potential electrodes. A) Average of the total accumulated charge on the positive side as a function of the applied potential. The error bars correspond to the standard deviation of the charge distribution. B) Normalized distributions of the instantaneous total charge on the positive half-electrode for various finite electric fields (E = 0.01, 0.1 and 1 V nm -1 from left to right, red lines) or the corresponding applied voltages (∆Ψ = 0.11, 1.1 and 11 V from left to right, blue lines). C) Snapshot of the positive electrode for a finite field of 0.5 V nm -1 (left) and the corresponding applied potential of 5.5 V (right), for the same electrolyte configuration. The atoms are colored according to their instantaneous charge, between 0 e (red) and 0.1 e (blue).

FIG. 4 :

 4 FIG.4: Charge density (A), electric field (B) and Poisson potential (C) profile for an applied field of 0.1V/nm (corresponding to a potential drop across the cell of 1.1 V. The gray-shaded area corresponds to the electrodes (the CCS was translated in order to match with the ECS setup).

  FIG. 7:Capacitance obtained using the fluctuationdissipation theorem for the finite electric field (red) and constant potential (blue). The grey-shaded area corresponds to the confidence interval obtained through the Q(∆Ψ) relation

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 771294).