δ¹⁸O of atmospheric oxygen (δ¹⁸O_atm) undergoes millennial-scale variations during the last glacial period, and systematically increases during Heinrich stadials (HSs). Changes in δ¹⁸O_atm combine variations in biospheric and water cycle processes. The identification of the main driver of the millennial variability in δ¹⁸O_atm is thus not straightforward. Here, we quantify the response of δ¹⁸O_atm to such millennial events using a freshwater hosing simulation performed under glacial boundary conditions. Our global approach takes into account the latest estimates of isotope fractionation factor for respiratory and photosynthetic processes and make use of atmospheric water isotope and vegetation changes. Our modeling approach allows to reproduce the main observed features of a HS in terms of climatic conditions, vegetation distribution and δ¹⁸O of precipitation. We use it to decipher the relative importance of the different processes behind the observed changes in δ¹⁸O_atm. The results highlight the dominant role of hydrology on δ¹⁸O_atm and confirm that δ¹⁸O_atm can be seen as a global integrator of hydrological changes over vegetated areas.