  We calculate the generalized soft functions at $\mathcal{O}(\alpha_s^2)$ at
next-to-leading power accuracy for the Drell-Yan process at threshold. The
operator definitions of these objects contain explicit insertions of soft gauge
and matter fields, giving rise to a dependence on additional convolution
variables with respect to the leading power result. These soft functions
constitute the last missing ingredient for the validation of the bare
factorization theorem to NNLO accuracy. We carry out the calculations by
reducing the soft squared amplitudes into a set of canonical master integrals
and we employ the method of differential equations to evaluate them. We retain
the exact $d$-dimensional dependence of the convolution variables at the
integration boundaries in order to regulate the fixed-order convolution
integrals. After combining our soft functions with the relevant collinear
functions, we perform checks of the results at the cross-section level against
the literature and expansion-by-regions calculations, at NNLO and partly at
N$^3$LO, finding agreement.
