We present the first calculation of the next-to-next-to-leading order threshold soft function for top quark pair production at hadron colliders, with full velocity dependence of the massive top quarks. Our results are fully analytic, and can be entirely written in terms of generalized polylogarithms. The scale-dependence of our result coincides with the well-known two-loop anomalous dimension matrix including the three-parton correlations, which at the two-loop order only appear when more than one massive partons are involved in the scattering process. In the boosted limit, our result exhibits the expected factorization property of mass logarithms, which leads to a consistent extraction of the soft fragmentation function. The next-to-next-to-leading order soft function obtained in this paper is an important ingredient for threshold resummation at the next-to-next-to-next-to-leading logarithmic accuracy.