In this paper, we investigate the solvability of a boundary value problem for a heat and mass transfer model with the spatially averaged Rayleigh function. The considered model describes the 3D steady-state non-isothermal flow of a generalized Newtonian fluid (with shear-dependent viscosity) in a bounded domain with Lipschitz boundary. The main novelty of our work is that we do not neglect the viscous dissipation effect in contrast to the classical Boussinesq approximation, and hence, deal with a system of strongly nonlinear partial differential equations. Using the properties of the averaging operation and d-monotone operators as well as the Leray–Schauder alternative for completely continuous mappings, we prove the existence of weak solutions without any smallness assumptions for model data. Moreover, it is shown that the set of all weak solutions is compact, and each solution from this set satisfies some energy equalities.