Classical MD simulation models were used to equilibrate the bulk solution and ice models of water. The bulk aqueous system consisted of 200 water molecules in a cubic simulation box (18.27 Å) with periodic boundary conditions (PBC) along all three dimensions. The equations of motion are integrated with the velocity Verlet algorithm with a time step of 0.5 fs. Simulations are carried out in the constant volume canonical ensemble (NVT), where the temperature is controlled by stochastic velocity rescaling with a relaxation time of 1 ps, and long-range electrostatics are computed using the particle-particle/particle-mesh solver (PPPM). MD simulations are performed using the LAMMPS package. Water molecules are modeled using the fixed-charge TIP4P/Ice force field, which was parameterized to reproduce the experimental melting temperature of ice at ambient conditions, and it accurately reproduces the equations of state of water and ice near the melting temperature and at mild supercooling. For the air-ice interface, a proton-disordered hexagonal ice model containing 192 water molecules was generated using a Monte Carlo procedure to minimize the total electric dipole, and it was cut into a slab of 20 Å thick along the basal plane. The orthorhombic simulation box measures 18.000 × 15.590 × 90.000 Å3 and periodic boundary conditions (PBC) are applied in all three directions. To avoid finite-size artifacts, periodic images of the slab are separated by ∼ 70 Å vacuum along the z direction. The system was equilibrated for 10 ns at 263 K to obtain a stable QLL at the surface. The QLL is about two ice bilayers thick at both surfaces, in agreement with previous experimental observations and MD simulations. First-Principles Molecular Dynamics FPMD simulations of several tens of picoseconds are performed for molecules in aqueous solution and at the air-ice interface. The purpose of these simulations is to attain a clearer picture of the dynamical hydrogen bonding network and to sample the configuration of the solvated molecule with a computational framework consistent with that used for the excited state calculations. DFT-based FPMD simulations were carried out using the Quickstep approach implemented in the CP2K 8.2 package, employing the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) for the exchange and correlation (XC) functional with D3 dispersion corrections. Valence Kohn-Sham orbitals are expanded on a double-ζ localized basis set 38 in real space, and core states are treated implicitly using Goedecker-Teter-Hutter pseudopotentials. Plane waves up to a cutoff energy of 300 Ry are used as a basis set for the density in reciprocal space. Using the aqueous and ice systems generated with classical MD and the aforementioned empirical force fields, the nitrate molecules and their counterions were placed at the surface and in the bulk. Eventually, using CP2K, the systems were equilibrated for 10 ps at the target temperature, followed by 70 ps long production runs in the constant volume canonical ensemble (NVT) where the temperature is controlled by stochastic velocity rescaling with a relaxation time of 1 ps. Aqueous solutions were run at room temperature (300 K) and ice slabs at 263 K. In these simulations, hydrogen atoms were replaced with deuterium, thus allowing a relatively large timestep of 0.5 fs to integrate the equations of motion. The structural properties reported in the results section, including density profiles, density maps, radial distribution functions (RDF), and coordination numbers, are calculated averaging over 1400 frames taken every 50 fs for each 70 ps trajectory. Uncertainties are computed as the standard deviation over seven 10-ps-long blocks into which each trajectory is split. 800 statistically independent molecular configurations were extracted from each trajectory to compute UV-visible absorption spectra with the ensemble method. The necessity to run long simulations, of the order of several hundreds of ps, for systems containing up to 192 water molecules limits our choice of XC functional to semilocal GGA. Several studies compared the accuracy of various XC functionals when modeling water, showing that standard GGA functionals lead to substantial overstructuring of the liquid and a large underestimate of its equilibrium density. However, the inclusion of dispersion corrections drastically improves the description of liquid water at nearly no additional computational cost. The recently developed meta-GGA strongly-constrained and appropriately-normed (SCAN)