Second Harmonic Generation (SHG) is a nonlinear optical process based on the conversion of two photons of the same frequency into a single photon of twice the fundamental frequency. Considering its high sensitivity to the inversion-symmetry breaking, SHG has emerged as one of the most powerful techniques used to selectively monitor surface dynamics and reactions for all types of interfaces as well as for imaging non-centrosymmetric structures. Nevertheless, the origin of the SHG signal at a molecular scale is poorly understood and its interpretation as well as its simulation still mainly rely on molecular hyperpolarizability and orientation. To gain understanding in the origin of SHG signal, we develop here a new ab initio method for the simulation of the experimental SHG response, based on the direct calculation of the surface second-order electric susceptibility that controls the SHG response. This method is freed from the hyperpolarizability concept, which strongly depends upon the local chemical environment of the molecules. A simple water-air interface is used as a case study to develop and validate the method. We demonstrate that the bulk quadrupole contribution significantly overwhelms the interface dipole term in the purely interfacial induced second-order polarization. By including the bulk quadrupole contribution in the calculation, we obtain simulated SHG responses that are in unprecedented agreement with the experimental signal, with average relative errors lower than 5%. This work highlights the modest influence of the molecular orientation and the high significance of the bulk quadrupole contribution, which, by definition, does not depend on the interface, in the total experimental response.