Reynolds Averaged Navier-Stokes (RANS) turbulence models have historically been developed using a combination of theoretical/physical/mathematical arguments, modeler expertise, and empirical data-fitting. As a consequence of the loss of information incurred during the ensemble averaging process, RANS model development has always been data-driven, to an extent, out of necessity. With the profusion of high resolution simulations and experimental methods, there exists a prodigious opportunity to more comprehensively inform closure models with data and improve their utlility in practical predictions of complex turbulent flows. In contrast to inferring model parameters, this work uses inverse modeling to directly obtain corrective, spatially distributed terms that are consistent with the predictive RANS models, thus offering a route to address the issue of model discrepancy which has long plagued turbulence model development. The posterior distributions of the inferred terms are quantified using Bayesian inversion and sampling. Results are presented for turbulent channel flow and flows with wall curvature. While the source terms inferred in this work are problem dependent, an ongoing subset of our work is to use machine learning methodologies to obtain general functional forms of the type that are inferred in this work. This paper specifically focuses on the mechanics of Bayesian inversion and its implications on the quantification of errors and uncertainties in turbulence modeling.