The interaction between an electric field and the electric charges in a material is described by electrostatic screening, which in metallic systems is commonly thought to be confined within a distance of the order of the Thomas-Fermi length. The validity of this picture, which holds for surface charges up to $sim 10^{13},mathrm{cm^{-2}}$, has been recently questioned by several experimental results when dealing with larger surface charges, such as those routinely achieved via the ionic gating technique. Whether these results can be accounted for in a purely electrostatic picture is still debated. In this work, we tackle this issue by calculating the spatial dependence of the charge carrier density in thin slabs of niobium nitride via an ab initio density functional theory approach in the field-effect transistor configuration. We find that perturbations induced by surface charges $lesssim 10^{14},mathrm{cm^{-2}}$ are mainly screened within the first layer, while those induced by larger surface charges $sim 10^{15},mathrm{cm^{-2}}$ can penetrate over multiple atomic layers, in reasonable agreement with the available experimental data. Furthermore, we show that a significant contribution to the screening of large fields is associated not only to the accumulation layer of the induced charge carriers at the surface, but also to the polarization of the pre-existing charge density of the undoped system.