Turbulent mixing in a methane–oxygen recessed injector is studied using direct numerical simulations. The operating point is chosen to be fuel-rich and at high pressure to recreate a representative environment for space propulsion applications. The results are used to investigate the transport of the turbulent mixture fraction statistics and the validity of conventional transport models. It is observed that molecular diffusion is only relevant near the boundary layer of the injection recess cavity and at the recirculation zone. Moreover, turbulent mixing in the axial direction is negligible as radial turbulent diffusion dominates. Radial turbulent diffusion near injection is driven by Kelvin–Helmholtz Instabilities (KHI) manifesting at large scales in the order of the injector geometry. The dominance of this process over microscale mixing originates negative turbulent diffusion, which produces a mixture resegregation and the appearance of lean pockets far from the oxidizer injection plane. Gradient models display poor capabilities for the prediction of this sort of phenomena. Closure models for the turbulent mixing transport terms are proposed and evaluated. An anisotropic gradient model is devised, providing performance improvements within the recess cavity and the recirculation region. In addition, a novel filtered Reynolds-averaged Navier Stokes approach based on the mixing state is proposed. This new methodology shows excellent prediction capabilities in the regions dominated by KHI, accurately predicting negative turbulent diffusivity. The challenges associated with this model are commented on, and strategies to enable its application are proposed.