We present a novel, data-driven closure model for Reynolds-averaged Navier-Stokes(RANS) equations which consists of two-parts. A parametric one, which a tensorbasis neural-network and a non-parametric one which makes use of latent, randomvariables in order to capture aleatoric model uncertainty. Our fully Bayesianformulation, incorporating sparsity-inducing priors, identifies areas of the problemdomain where the parametric closure falls short, requiring stochastic corrections tothe Reynolds stress tensor. Training employs sparse, indirect data such as meanvelocities and pressures, in contrast to the majority of alternatives which requiredirect, Reynolds stress data. For inference and learning, we employ StochasticVariational Inference, facilitated by an adjoint-based differentiable solver. Thisend-to-end differentiable framework can ultimately yield accurate, probabilisticpredictions for flow quantities, even in regions with model errors, as exemplifiedby the backward-facing step benchmark problem.