Ocean general circulation models (OGCMs) contain numerous parameterizations of sub-grid scale processes. The parameter tuning procedure is rarely reported and often done by hand. We present an automated alternative: Bayesian optimization, a method which has recently emerged as a frontier in expensive black box optimization. VerOpt, a Python package for the ocean model Veros, adapts Bayesian optimization to climate model tuning. We use VerOpt to identify a set of parameter values of the Turbulent Kinetic Energy (TKE) closure scheme that minimize mixed layer depth (MLD) bias in Veros. We present the results of two optimization procedures: TWIN and OBS. The goal is to minimize modeled MLD error relative to a target map. In TWIN, the target is MLD simulated using Veros with a known parameterization. The ratio of two TKE parameters ckcϵ-1, proportional to the critical Richardson number Ric, is the dominant factor in setting the global MLD. After 180 model simulations, the lowest error in the TWIN experiment is 1.18%. In OBS, the target is MLD climatology. The MLD bias is smallest when Ric<1, and the default TKE parameterization falls within this range. We find, however, that altering the TKE parameterization is not sufficient to reduce the significant MLD bias of 42.62%. The OBS experiment results indicate that the TKE scheme parameters are not the dominant source of MLD bias in Veros. We discuss other possible sources of MLD bias, as well as the potential of extending of the optimization procedure to other parameterizations.