This paper proposes a density-based topology optimization method for natural convection problems using the lattice Boltzmann method (LBM). As the LBM can be developed as a completely explicit scheme, its attractive features over the traditional ones, such as the finite element method, are (1) suitability for solving unsteady flow problems and (2) scalability for large-scale parallel computing. We develop an LBM code for solving unsteady natural convection problems and provide its sensitivity analysis based on the so-called adjoint lattice Boltzmann method. Notably, the adjoint equation is derived from the discrete particle velocity Boltzmann equation and can be solved similarly to the original LBM concerning unsteady natural convection problems. We first show that the proposed method can produce similar results to the previous work in a steady-state natural convection problem. We then demonstrate the efficacy of the proposed method through 2D numerical examples concerning unsteady natural convection. As a large-scale problem, we tackle a 3D unsteady natural convection problem on a parallel supercomputer. All the developed codes written in C++ are available at https://github.com/PANFACTORY/PANSLBM2.git.