Abstract. A multiscale modeling approach for studying the ocean surface turbulent mixing is explored by coupling an ocean general circulation model (GCM) MPAS-Ocean with the PArallel Large eddy simulation Model (PALM). The coupling approach is similar to the superparameterization approach that has been used mostly to represent the effects of deep convection in atmospheric GCMs. However, since the emphasis here is on the small-scale turbulent mixing processes and their interactions with the larger-scale processes, a high-fidelity, three-dimensional large eddy simulation (LES) model is used, in contrary to a simplified process-resolving model with reduced physics or reduced dimension commonly used in the superparameterization literature. To reduce the computational cost, a customized version of PALM is ported on the general-purpose graphics processing unit (GPU) with OpenACC, achieving 10–16 times overall speedup as compared to running on a single CPU. Even with the GPU-acceleration technique, superparameterization of the ocean surface turbulent mixing using high-fidelity and three-dimensional LES over the global ocean in GCMs is still computationally intensive and infeasible for long simulations. However, running PALM regionally on selected MPAS-Ocean grid cells is shown to be a promising approach moving forward. The flexible coupling between MPAS-Ocean and PALM outlined here allows further exploration of the interactions between ocean surface turbulent mixing and larger-scale processes, and development of better ocean surface turbulent mixing parameterizations in GCMs.