An Euler–Lagrange multicomponent, non-Newtonian Lattice-Boltzmann method is applied for the first time to model a full-scale gas-mixed anaerobic digester for wastewater treatment. Rheology is modelled through a power-law model and, for the first time in gas-mixed anaerobic digestion modelling, turbulence is modelled through a Smagorinsky Large Eddy Simulation model. The hydrodynamics of the digester is studied by analysing flow and viscosity patterns, and assessing the degree of mixing through the Uniformity Index method. Results show independence from the grid size and the number of Lagrangian substeps employed for the Lagrangian sub-grid simulation model. Flow patterns are shown to depend mildly on the choice of bubble size, but not the asymptotic degree of mixing. Numerical runs of the model are compared to previous results in the literature, from a second-ordered Finite-Volume Method approach, and demonstrate an improvement, compared to literature data, of 1000-fold computational efficiency, massive parallelizability and much finer attainable spatial resolution. Whilst previous research concluded that the application of LES to full-scale anaerobic digestion mixing is unfeasible because of high computational expense, the increase in computational efficiency demonstrated here, now makes LES a feasible option to industries and consultancies.