<p>Numerical modeling of the process of combined heat transfer by natural convection and thermal surface radiation was performed. The velocity vector components were calculated from the mesoscopic distribution function of conditional particles within the framework of the lattice Boltzmann method. On the other hand, the calculation of the thermodynamic characteristics of the thermogravitational flow was carried out by a finite-difference solution of the macroscopic energy equation. The developed hybrid algorithm, simulating the direct numerical modeling method, was tested on reference experimental data on turbulent natural convection in a closed rectangular cavity that were obtained by other researchers. The range of variation of the Rayleigh number, conductive-radiative parameter and wall emissivity was 4⋅10⁴ ≤ Ra ≤ 10⁹, 55.84 ≤ <i>N</i>_r ≤ 163.448 and 0 ≤ ε ≤ 1, respectively. In the course of numerical modeling it was established that the process of convective-radiative heat transfer in a closed square cavity with lateral heating/cooling becomes nonstationary at Ra ≥ 6⋅10⁷ and <i>N</i>_r ≥ 63.95. In this case, the evolution trend of the average convective and effective Nusselt numbers acquires an oscillatory character with a constant amplitude of oscillations in the range of 6⋅10⁸ ≤ Ra ≤ 10⁹, 63.95 ≤ <i>N</i>_r ≤ 75.89.</p>