A bstract We consider QCD radiative corrections to Higgs boson pair production through gluon fusion in proton collisions. We combine the exact next-to-leading order (NLO) contribution, which features two-loop virtual amplitudes with the full dependence on the top quark mass M t , with the next-to-next-to-leading order (NNLO) corrections computed in the large- M t approximation. The latter are improved with different reweighting techniques in order to account for finite- M t effects beyond NLO. Our reference NNLO result is obtained by combining one-loop double-real corrections with full M t dependence with suitably reweighted real-virtual and double-virtual contributions evaluated in the large- M t approximation. We present predictions for inclusive cross sections in pp collisions at s = 13, 14, 27 and 100 TeV and we discuss their uncertainties due to missing M t effects. Our approximated NNLO corrections increase the NLO result by an amount ranging from +12% at s = 13 TeV to +7% at s = 100 TeV, and the residual uncertainty of the inclusive cross section from missing M t effects is estimated to be at the few percent level. Our calculation is fully differential in the Higgs boson pair and the associated jet activity: we also present predictions for various differential distributions at s = 14 and 100 TeV, and discuss the size of the missing M t effects, which can be larger, especially in the tails of certain observables. Our results represent the most advanced perturbative prediction available to date for this process.