Here we employ high-level quantum chemistry methods to compute core-level IE of bulk water. In addition to providing the best theoretical estimate of the bulk IE, we also aim to carefully investigate the convergence of the spectrum with respect to the details of computational protocol, to aid future theoretical studies. We use MD simulations with classical forcefields and ab initio potentials to simulate bulk water, and then use the snapshots from the MD simulations to carry out QM/MM (quantum mechanics/molecular mechanics) cal- with triple excitations). The main challenge in applying these methods to modeling condensed-phase phenomena is how to properly account for the effect of the solvent via embedding. Here we show that in calculations of core-level IEs simple electrostatic embedding (QM/MM) 47,48 converges slowly with respect to the size of the QM system, which illustrates the high sensitivity of the core-level states to the environment. We also investigate the contribution of different types of structures present in bulk water in the overall spectrum.