In this paper we consider the problem of model-based learning of turbulent flows using mobile robots. The key idea is to use empirical data to improve on numerical estimates of time-averaged flow properties that can be obtained using Reynolds-Averaged Navier Stokes (RANS) models. RANS models are computationally efficient and provide global knowledge of the flow but they also rely on simplifying assumptions and require experimental validation. In this paper, we instead construct statistical models of the flow properties using Gaussian Processes (GPs) and rely on the numerical solutions obtained from RANS models to inform their mean. We then utilize Bayesian inference to incorporate empirical measurements of the flow into these GPs, specifically, measurements of the time-averaged velocity and turbulent intensity fields. Our formulation accounts for model ambiguity and parameter uncertainty via hierarchical model selection. Moreover, it accounts for measurement noise by systematically incorporating it in the GP models. To obtain the velocity and turbulent intensity measurements, we design a cost-effective mobile robot sensor that collects and analyzes instantaneous velocity readings. We control this mobile robot through a sequence of waypoints that maximize the information content of the corresponding measurements. The end result is a posterior distribution of the flow field that better approximates the real flow and also quantifies the uncertainty in the flow properties. We present experimental results that demonstrate considerable improvement in the prediction of the flow properties compared to pure numerical simulations.