The lattice Bhatnagar-Gross-Krook (LBGK) model has become the most popular one in the lattice Boltzmann method for simulating the convection heat transfer in porous media. However, the LBGK model generally suffers from numerical instability at low fluid viscosities and effective thermal diffusivities. In this paper, a modified LBGK model is developed for incompressible thermal flows in porous media at the representative elementary volume scale, in which the shear rate and temperature gradient are incorporated into the equilibrium distribution functions. With two additional parameters, the relaxation times in the collision process can be fixed at a proper value invariable to the viscosity and the effective thermal diffusivity. In addition, by constructing a modified equilibrium distribution function and a source term in the evolution equation of temperature field, the present model can recover the macroscopic equations correctly through the Chapman-Enskog analysis, which is another key point different from previous LBGK models.Several benchmark problems are simulated to validate the present model with the proposed local computing scheme for the shear rate and temperature gradient, and the numerical results agree well with analytical solutions and/or those well-documented data in previous studies. It is also shown that the present model and the computational schemes for the gradient operators have a second-order accuracy in space, and better numerical stability of the present modified LBGK model than previous LBGK models is demonstrated. and mineral processing. Over the past several decades, considerable investigations and applications have been devoted to the convection heat transfer in porous media by many researchers using various traditional numerical methods (eg. finite difference, element and volume methods). A comprehensive review on this subject has been made by Cheng [1], Nield and Bejan [2] and Vafai [3].As a powerful computational tool based on the kinetic theory, the lattice Boltzmann method (LBM) has been successfully applied to simulating complex fluid flows and modeling the physics in fluids [4,5,6]. Some attractive advantages of LBM over the traditional numerical methods have been explained in references [7] and [8]. Thanks to the mesoscopic and kinetic nature, the LBM has been widely applied to the fluid flow and thermal problems in porous media [9, 10, 11] after its emergence [12]. Generally speaking, applications of the LBM to porous flows in the literature center around two scales: the pore scale and the representative elementary volume (REV) scale. At the pore scale [12,13,14,15,16,17], the standard lattice Boltzmann equation (LBE) is used to simulate fluid flows in pores, and the local information of flow can be directly obtained. Therefore, the LBM at the pore scale can be severed as the most straightforward way to investigate the macroscopic relations and reveal the microscopic mechanism of porous flows. However, the detailed geometric information of the pores is needed for this approach,...