While various parameterizations of vertical turbulent fluxes at different levels of complexity have been proposed, each has its own limitations. For example, simple first-order closure schemes such as the K-Profile Parameterization (KPP) lack energetic constraints; two-equation models like $k-\epsilon$ directly solve an equation for the turbulent kinetic energy but do not account for non-local fluxes, and high-order closures that include the non-local transport terms are computationally expensive. To address these, here we extend the Assumed-Distribution Higher-Order Closure (ADC) framework originally proposed for the atmospheric boundary layer and apply it to the ocean surface boundary layer (OSBL). By assuming a probability distribution function relationship between the vertical velocity and tracers, all second-order and higher-order moments are exactly constructed and turbulence closure is achieved in the ADC scheme. In addition, the ADC parameterization scheme has full energetic constraints. We have tested the ADC scheme against a combination of large eddy simulation (LES), KPP, and $k-\epsilon$ for surface buoyancy-driven convective mixing and found that the ADC scheme is robust with different vertical resolutions and compares well to the LES results.