Inaccuracies in the models of the Earth system, i.e., structural and parametric model errors, lead to inaccurate climate change projections.Errors in the model can originate from unresolved phenomena due to a low numerical resolution, as well as misrepresentations of physical phenomena or boundaries (e.g., orography). Therefore, such models lead to inaccurate short–term forecasts of weather and extreme events, and more importantly, long term climate projections. While calibration methods have been introduced to address for parametric uncertainties, e.g., by better estimation of system parameters from observations, addressing structural uncertainties, especially in an interpretable manner, remains a major challenge.Therefore, with increases in both the amount and frequency of observations of the Earth system, algorithmic innovations are required to identify interpretable representations of the model errors from observations. We introduce a flexible, general-purpose framework to discover interpretable model errors, and show its performance on a canonical prototype of geophysical turbulence, the two–level quasi–geostrophic system. Accordingly, a Bayesian sparsity–promoting regression framework is proposed, that uses a library of kernels for discovery of model errors. As calculating the library from noisy and sparse data (e.g., from observations) using convectional techniques leads to interpolation errors, here we use a coordinate-based multi–layer embedding to impute the sparse observations. We demonstrate the importance of alleviating spectral bias, and propose a random Fourier feature layer to reduce it in the proposed embeddings, and subsequently enable an accurate discovery. Our framework is demonstrated to successfully identify structural model errors due to linear and nonlinear processes (e.g., radiation, surface friction, advection), as well as misrepresented orography.