Hydration repulsion dominates the interaction between polar surfaces in water at nanometer separations and ultimately prevents the sticking together of biological matter. Although confirmed by a multitude of experimental methods for various systems, its mechanism remained unclear. A simulation technique is introduced that yields accurate pressures between solvated surfaces at prescribed water chemical potential and is applied to a stack of phospholipid bilayers. Experimental pressure data are quantitatively reproduced and the simulations unveil a rich microscopic picture: Direct membrane-membrane interactions are attractive but overwhelmed by repulsive indirect water contributions. Below about 17 water molecules per lipid, this indirect repulsion is of an energetic nature and due to desorption of hydration water; for larger hydration it is entropic and suggested to involve water depolarization. This antagonistic nature and the presence of various compensating contributions indicate that the hydration repulsion is less universal than previously assumed and rather involves finely tuned surface-water interactions.solvation | MD simulation | phospholipids H ydration repulsion (HR) universally acts between wellsolvated surfaces in water and balances the van der Waals attraction in the nanometer range. It ultimately prevents the collapse of biological matter and thereby provides macromolecular assemblies with the necessary lubrication for vital functioning, even in the congested cell environment. Although complex in its nature, it is rightfully considered a fundamental force in solution chemistry and structural biology (1). HR was first quantified experimentally for stacks of charge-neutral phospholipid bilayer membranes in terms of pressure-distance curves (2-4), confirmed for two individual bilayers by the surface force apparatus (SFA) (5, 6), and is now known to universally act between nucleic acids, proteins, and polysaccharides alike (7). It exhibits an exponential decay with a decay length of a few Ångstrom (4) and as a heuristic law is nowadays commonly used in modeling the forces between polar surfaces in water (8). Although several theoretical (9-11) and simulation (12-18) studies elucidated partial aspects of the HR, none treated the full complexity of the problem and could quantitatively reproduce and explain experimental pressure-distance curves, meaning that the HR mechanism remained essentially unclear. The reason for this is obvious: Theory typically only treats one part of the problem, be it the water-water interactions, the water-surface binding, or the configurational entropy of bilayer molecules, whereas current simulation strategies account for the constant water chemical potential either in the form of a large reservoir (13-15) or by grand-canonical simulations (16,17). Due to limitations in the numerical accuracy, however, both approaches do not enable quantitative comparison of the HR pressure with experimental data. We solve this problem by introducing the thermodynamic extrapolation method (TE...