In this talk the projection-based reduced order model (ROMs) in a context of variational multiscale (VMS) formulation for modeling wind turbines will be presented. The flow is modeled using the Navier-Stokes equations for incompressible stratified flows with the Boussinesq approximation, incorporating an additional advection-diffusion equation for the temperature/density field. Finite-element method is used for spatial discretization of the problem and VMS formulation is adopted for turbulence modeling. For ROM, the POD-Galerkin method with mesh-based hyper reduction is used. The wind turbine is modeled using the actuator line model (ALM). The accuracy and robustness of the formulation is shown on various 3D cases including NREL 5MW wind turbine, NREL Phase VI wind turbine and simulation of multiple NREL 5MW wind turbines. For all applications the Reynolds number is an order of millions. Comparisons with field data, laboratory experiments, and high-fidelity model demonstrate the robustness of the consistent ROM-ALM-VMS formulation.