In uncertainty quantification of computational models (e.g., turbulence modeling) with Bayesian inferences, Gaussian processes are commonly used as the prior for model discrepancies. However, constructing the covariance kernel is a challenging task that requires significant physical knowledge of the problem. On the other hand, often the model discrepancies are described by partial differential equations (PDEs) of known structures (e.g. Reynolds stress transport equations for turbulent flows), albeit with unclosed terms (e.g, velocity triple correlation and press–strain-rate correlation). In this work, we utilize such PDEs to construct physics-informed covariance kernels by exploiting the fundamental connection between PDEs and covariance functions. We demonstrate the merits of the physics-informed covariance kernel with the application to turbulent flows over periodic hills. The covariance kernel and the modes (eigenfunctions) obtained from the physics-informed approach are physically more realistic than those obtained with commonly used squared exponential kernels. This method also has the potential of improving the performance of Bayesian model-form uncertainty quantification in applications beyond turbulent flows.