Molecular dynamics simulations were performed using the LAMMPS interface of the n2p2 code. Simulation cells were periodic in all three directions. In order to avoid spurious interactions between periodic images of the slab in the z-direction (perpendicular to the graphene sheets), we included a layer of vacuum at least three times the slit height between periodic images of the slab and applied the Yeh and Berkowitz correction. Coulombic interactions were evaluated with the PPPM method. Equations of motion were integrated using the velocity-Verlet algorithm with a time step of 0.5 fs. The interatomic forces used for this numerical integration are based on the average prediction of all eight committee members of the neural network potential. All simulations were carried out in the NVT ensemble using a Nosé-Hoover style thermostat with a damping parameter of 100 fs. We chose to run all simulations at a temperature of 350K to avoid the complicating influence of phase transitions in the bilayer electrolyte, which were observed at 300K (see the SI for more details). Simulation cells were initially prepared by randomly packing the desired number of water molecules and ions between two graphene sheets. These initial structures were subsequently equilibrated with a 1 ns simulation using a classical force field (SPC/E water model, NaCl parameters from Dang, 88 and water-carbon parameters from Werder ), followed by a 1 ns simulation using the NNP. Production runs were then carried out for 5 ns. For each system, we performed five independent replicate simulations. All error bars in the text reflect the standard deviation among these replicates.