Compressed Natural Gas (CNG) is a promising alternative fuel for internal combustion engines as its combustion is fuel-efficient and lean in carbon dioxide compared to gasoline. The high octane number of methane gives rise to significant increase of the thermodynamic efficiency due to higher possible compression ratios. In order to use this potential, new stratified mixture formation concepts for CNG are investigated by means of numerical fluid simulations. For decades RANS methods have been the industry standard to model three-dimensional flows. Indeed, there are well-known deficiencies of the widely used eddy viscosity turbulence models based on the applied Boussinesq hypothesis. Reynolds stress turbulence models as well as scale resolving simulation approaches can be appealing alternative choices since they offer higher accuracy. However, due to their large computing effort, they are still mostly impractical for the daily use in industrial product development processes. A more suitable solution seems to be the application of an explicit algebraic Reynolds stress (EARSM) turbulence model based on a non-linear relation between the Reynolds stresses and the mean strain-rate. This allows a better prediction of various flow effects in the simulations. In the presented work the direct injection of CNG has been simulated by means of commercial RANS-based CFD using various turbulence models. The simulation set-up includes the discretization of one cylinder and part of the centrally mounted injector geometry. High pressure ratios are present in between the rail and the cylinder resulting in a supersonic gas jet emerging from the injector and mixing into the surrounding air. Mass transfer in between both components is modeled by a passive scalar transport equation based on the gradient diffusion hypothesis using a constant global turbulent Schmidt number. It is shown in this work that the capturing of anisotropic turbulence effects is crucial for the correct modeling of the mixture formation because turbulent diffusion is dominant in the mixing layer of supersonic gas jets. Overall an explicit algebraic Reynolds stress model predicts more pronounced mixing and shows better agreement with experimental data than the commonly used eddy viscosity turbulence models.