Bu↵ered Verlet lists were used to keep track of atomic neighbours. Electrostatic interactions were calculated using the smooth Particle Mesh Ewald algorithm. For the Lennard-Jones potential, we used di↵erent cuto↵ values of 0.8, 1.0, 1.2, and 1.4 nm. CHARMM force fields use a switch function for the LJ potential, but this would introduce an extra parameter-the distance at which the switching begins-and thus we decided to always shift the potential to zero at the cuto↵. We applied dispersion corrections to energy and pressure, as these correction is used for monolayer simulations with CHARMM. However, the e↵ect of dispersion corrections on the water-air surface tension is within the error estimate. We also repeated the simulations using LJ-PME. In all simulations, temperature was controlled by the stochastic velocity rescaling algorithm with a target temperature of either 298, 310, or 323 K, and a time constant of 1 ps. The geometry of the water molecules was constrained by the SETTLE algorithm. The surface tension values were extracted from pressure components normal (P N ) and lateral (P L ) to the interface as