We present the two-loop virtual corrections to gg→ZH involving a closed top-quark loop. These amplitudes provide the missing ingredient required to compute the NLO QCD corrections to Higgs boson production in association with a Z-boson in gluon fusion. The two-loop integrals, involving massive top-quarks, are first expressed in terms of finite integrals and are then calculated numerically. We will describe the techniques used to obtain and evaluate these amplitudes and present comparisons to various approximations in the literature. We will discuss an upcoming public release of pySecDec, the tool used to compute these corrections, and highlight a new feature that simplifies the numerical computation of amplitudes.