Abstract: Higgs-pair production via gluon fusion is the dominant production mechanism of Higgs-boson pairs at hadron colliders. In this work, we present details of our numerical determination of the full next-to-leading-order (NLO) QCD corrections to the leading top-quark loops. Since gluon fusion is a loop-induced process at leading order, the NLO calculation requires the calculation of massive two-loop diagrams with up to four different mass/energy scales involved. With the current methods, this can only be done numerically, if no approximations are used. We discuss the setup and details of our numerical integration. This will be followed by a phenomenological analysis of the NLO corrections and their impact on the total cross section and the invariant Higgs-pair mass distribution. The last part of our work will be devoted to the determination of the residual theoretical uncertainties with special emphasis on the uncertainties originating from the scheme and scale dependence of the (virtual) top mass. The impact of the trilinear Higgs-coupling variation on the total cross section will be discussed.