We compute the top quark contribution to the two-loop amplitude for on-shell $Z$ boson pair production in gluon fusion, $gg to ZZ$. Exact dependence on the top quark mass is retained. For each phase space point the integral reduction is performed numerically and the master integrals are evaluated using the auxiliary mass flow method, allowing fast computation of the amplitude with very high precision.