  We make a comparison of two tadpole renormalization schemes in the context of
the quarkonium hyperfine splittings in lattice NRQCD. Improved gauge-field and
NRQCD actions are analyzed using the mean-link $u_{0,L}$ in Landau gauge, and
using the fourth root of the average plaquette $u_{0,P}$. Simulations are done
for $c\bar c$, $b\bar c$, and $b\bar b$ systems. The hyperfine splittings are
computed both at leading and at next-to-leading order in the relativistic
expansion. Results are obtained at lattice spacings in the range of about
0.14~fm to 0.38~fm. A number of features emerge, all of which favor tadpole
renormalization using $u_{0,L}$. This includes much better scaling behavior of
the hyperfine splittings in the three quarkonium systems when $u_{0,L}$ is
used. We also find that relativistic corrections to the spin splittings are
smaller when $u_{0,L}$ is used, particularly for the $c\bar c$ and $b\bar c$
systems. We also see signs of a breakdown in the NRQCD expansion when the bare
quark mass falls below about one in lattice units. Simulations with $u_{0,L}$
also appear to be better behaved in this context: the bare quark masses turn
out to be larger when $u_{0,L}$ is used, compared to when $u_{0,P}$ is used on
lattices with comparable spacings. These results also demonstrate the need to
go beyond tree-level tadpole improvement for precision simulations.
