This study utilizes Deep Neural Networks (DNN) to improve the K-Profile Parameterization (KPP) for the vertical mixing effects in the ocean’s surface boundary layer turbulence. The DNNs were trained using 11-year turbulence-resolving solutions, obtained by running a large eddy simulation model for Ocean Station Papa, to predict the turbulence velocity scale coefficient and unresolved shear coefficient in the KPP. The DNN-augmented KPP schemes (KPP_DNN) have been implemented in the General Ocean Turbulence Model (GOTM). This implementation is stable for long-term integration and as efficient as existing variants of KPP schemes. Three different KPP_DNN schemes, varying in input and output variables, have been developed and trained. The performance of models using the KPP_DNN schemes is compared with that of those using popular deterministic first-order and second-moment closure turbulent mixing parameterizations. Solution comparisons show that the simulated mixed layer is cooler and deeper, aligning closely with observations when wave effects are included in parameterizations. In the KPP framework, changes to the velocity scale of unresolved shear, which is used to calculate mixed layer depth, have a larger impact on the simulated mixed layer than do changes to the magnitude of diffusivity. In the KPP_DNN, changes to unresolved shear depend on not only on wave forcing, but also on the mixed layer depth and buoyancy forcing.