We compute the topological charge and its susceptibility in finite temperature (2+1)-flavor QCD on the lattice applying a gradient flow method. With the Iwasaki gauge action and nonperturbatively $O(a)$-improved Wilson quarks, we perform simulations on a fine lattice with~$asimeq0.07,mathrm{fm}$ at a heavy $u$, $d$ quark mass with $m_pi/m_rhosimeq0.63$ but approximately physical $s$ quark mass with $m_{eta_{ss}}/m_phisimeq0.74$. In a temperature range from~$Tsimeq174,mathrm{MeV}$ ($N_t=16$) to $697,mathrm{MeV}$ ($N_t=4$), we study two topics on the topological susceptibility. One is a comparison of gluonic and fermionic definitions of the topological susceptibility. Because the two definitions are related by chiral Ward-Takahashi identities, their equivalence is not trivial for lattice quarks which violate the chiral symmetry explicitly at finite lattice spacings. The gradient flow method enables us to compute them without being bothered by the chiral violation. We find a good agreement between the two definitions with Wilson quarks. The other is a comparison with a prediction of the dilute instanton gas approximation, which is relevant in a study of axions as a candidate of the dark matter in the evolution of the Universe. We find that the topological susceptibility shows a decrease in $T$ which is consistent with the predicted $chi_mathrm{t}(T) propto (T/T_{rm pc})^{-8}$ for three-flavor QCD even at low temperature $T_{rm pc} < Tle1.5 T_{rm pc}$.