  We present a novel data-driven method for determining the hadronic
interaction strengths of axion-like particles (ALPs) with QCD-scale masses.
Using our method, it is possible to calculate the hadronic production and decay
rates of ALPs, along with many of the largest ALP decay rate to exclusive final
states. To illustrate the impact on QCD-scale ALP phenomenology, we consider
the scenario where the ALP-gluon coupling is dominant over the ALP coupling to
photons, electroweak bosons, and all fermions for $m_{\pi} \lesssim m_a
\lesssim 3$ GeV. We emphasize, however, that our method can easily be
generalized to any set of ALP couplings to SM particles. Finally, using the
approach developed here, we provide calculations for the branching fractions of
$\eta_c \to VV$ decays, i.e. $\eta_c$ decays into two vector mesons, which are
consistent with the known experimental values.
