Temperature dependence of the $U_A(1)$ anomaly is investigated by taking into account mesonic fluctuations in the $U(3)\times U(3)$ linear sigma model. A field dependent anomaly coefficient function of the effective potential is calculated within the finite temperature functional renormalization group approach. The applied approximation scheme is a generalization of the chiral invariant expansion technique developed in [G. Fejos, Phys. Rev. D 90, 096011 (2014)]. We provide an analytic expression and also numerical evidence that depending on the relationship between the two quartic couplings, mesonic fluctuations can either strengthen of weaken the anomaly as a function of the temperature. The role of the six-point invariant of the $U(3)\times U(3)$ group, and therefore the stability of the chiral expansion is also discussed in detail.