Parameterizations of unresolved turbulent processes in the ocean compromise the fidelity of large-scale ocean models used in climate change projections. In this work, we use a Bayesian approach for evaluating and developing turbulence parameterizations by comparing parameterized models with observations or high-fidelity numerical simulations. The method obtains optimal parameter values, correlations, sensitivities, and, more generally, likely distributions of uncertain parameters. We demonstrate the approach by estimating the uncertainty of parameters in the popular ‘K-Profile Parameterization’, using an ensemble of large eddy simulations of turbulent penetrative convection in the ocean surface boundary layer. We uncover structural deficiencies and discuss their cause. We conclude by discussing the applicability of the approach to Earth system models.