Deep learning has achieved noteworthy success in various applications. There is a new trend of adopting deep learning methods to study partial differential equations (PDEs). In this article, we use a kind of parallel deep neural networks (PDNNs) to approximate the solution of the Navier-Stokes equations coupled with the heat equation. The approach embeds the PDE formula into the loss function of the neural networks and the resulting networks is trained to meet the equations along with the boundary conditions. Moreover, the approximate ability of the PDNNs is demonstrated by the theoretical analysis of convergence.Finally, we present numerous numerical examples to verify effectiveness of the PDNNs, in which the viscosity 𝜈 is a constant, a function dependent on the space variable and a function of the temperature 𝜃.