Computational Fluid Dynamics based on RANS models remain the standard but suffer from high errors in complex flows. In particular, turbulent kinetic energy is over-produced in high strain rate regions, such as the near wake of wind turbine flows. Data-driven turbulence modelling methods aim to derive novel turbulence models with lower uncertainties, which generalize well to a certain class of flows. These state-of-the-art constitutes to first derive model-form corrections of a selected baseline model from high fidelity reference data, followed by regressing the corrections in terms of RANS-known flow features. For data-driven wind turbine wake modelling, industrial-scale wind turbines and non-neutral atmospheric boundary layers have yet to be considered. In this thesis, the first steps are made to address this research gap. First, Large-Eddy Simulation data is generated and validated against literature. The considered cases are under neutral, convective and stable atmospheric conditions. The frozen-RANS methodology, a technique used to derive turbulence model corrections given the high fidelity data, is then extended to for non-neutral conditions. The new framework now provides corrections to both the Boussinesq eddy viscosity hypothesis for the Reynolds stress and the gradient-diffusion hypothesis for the turbulent heat flux. By injecting the obtained corrections into dynamic RANS simulation, the baseline turbulence model deficiencies are corrected. In particular, high rate-of-strain regions now no longer show an overproduction of mechanical turbulence. Similarly, the lack of buoyant turbulence production in the free-stream atmosphere under convective conditions is solved. In the stable case, buoyant destruction is too large in the free-stream but not large enough in the wake. For the neutral and stable case, the corrected models produce wake velocity profiles that show excellent agreement with the large-eddy simulation reference data. Issues in the wall stress solution of the convective large-eddy simulation propagate to issues in the corrected RANS solutions, proving the necessity of high-quality data. Furthermore, it is shown that for most cases a single scalar correction to the turbulent heat flux, as opposed to the full vector correction, is sufficient for improving the error introduced by the gradient-diffusion hypothesis. This result is considerable since the simpler correction would be much easier to regress in terms of mean RANS-known quantities. The computational cost of the corrected RANS models is around the same as that of baseline RANS models; only 2\%-5\% of the large-eddy simulation computational cost.