This dissertation presents a study of hyperbolic turbulence modelling for the Gaussian ten-moment equations. In gaskinetic theory, moment closures offer the possibility of deriving a series of gas-dynamic governing equations from the Boltzmann equation. One typical example, the Gaussian ten-moment model, which takes the form of hyperbolic-relaxation equations, is considered as a competitive model for viscous gas flow when heat transfer effects are negligible. The hyperbolic nature of this model gives it several numerical advantages, compared to the Navier-Stokes equations. However, until this study, the application of the ten-moment equations has been limited to laminar flows, due to the lack of appropriate turbulence models. In this work, the ten-moment equations are, for the first time, Reynolds-averaged. The resulting equations inherit the hyperbolic balance-law form from the original equations with new unknowns, which require approximation by turbulence models. Most of the traditional turbulence models for the Reynolds-averaged Navier-Stokes equations are not perfectly well-suited for the Reynolds-averaged ten-moment equations, because the second-order derivatives presented in these models can break the pure hyperbolic nature of the original model. The relaxation methods are therefore proposed in this project to reform the existing turbulence models. Two relaxation methods, the Chen-Levermore-Liu p-system and Cattaneo-Vernotte models, are used to hyperbolize the Prandtl’s one-equation model, standard k-ε model and Wilcox k-ω model. The hyperbolic versions of these turbulence models are first shown to be equivalent to their original forms. They are then coupled to the Reynolds-averaged ten-moment equations to build the overall hyperbolic governing equations for turbulence flows. An axisymmetric version of Reynolds-averaged ten-moment equations is also derived. A dispersion analysis is conducted for the resulting governing equations, which shows the corresponding dispersive behaviour and stability. The effect of the relaxation parameters is investigated through several numerical tests. All derived turbulence models are applied to solve canonical validation test problems, including two-dimensional planar mixing-layer, free-jet and circular free-jet. The numerical evaluations are analysed and compared against existing experimental measurements.