To address this issue, FAU University together with IFPEN are currently implementing actuator-line based wind turbine models in the waLBerla LBM flow solver, with particular attention devoted to the collision operators. Whereas the Lattice-Boltzmann method is well-suited for lower Reynolds number regimes, standard collision operators suffer from instabilities at the high Reynolds numbers observed in wind energy applications. To obtain a stable and meaningful simulation, several investigations are conducted, typically using the Single and/or Multiple Relaxation Time collision operators together with different subgrid-scale models, such as e.g. the Smagorinsky and the WALE models. As this approach failed to stabilize the simulation sufficiently, more advanced collision operators, namely entropic, regularised, and cumulant-based are studied. For the wind turbine application, the cumulant Lattice-Boltzmann methods performed best in terms of stability. They have already been validated considering single or tandem wind turbines, and are being ported to Graphical Processing Units (GPU). The following figures illustrate the obtained results using CPU and GPU resources for the NewMexico wind Turbine. For this compressible simulation, Smagorinsky model is used with D3Q27 stencil and SRT collision operator. The simulation on one GPU (NVIDIA P100) substantially and expectedly outperforms the simulation on one CPU core (Intel Skylake G-6140) without any accuracy loss. Eight full nodes, i.e. 144 CPU cores, are required to obtain the same time-to-solution for both CPU and GPU systems.