  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~$a\simeq0.07\,\mathrm{fm}$ at a heavy $u$, $d$ quark mass with
$m_\pi/m_\rho\simeq0.63$ but approximately physical $s$ quark mass with
$m_{\eta_{ss}}/m_\phi\simeq0.74$. In a temperature range
from~$T\simeq174\,\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} < T\le1.5 T_{\rm pc}$.
