We present a systematic framework for the calculation of soft functions that are defined in
terms of $N\geq2$ light-like Wilson lines. The formalism represents an extension of a method
that we developed earlier for the calculation of dijet soft functions to the general $N$-jet
case. We discuss the technical aspects of this generalisation, focussing on SCET-1 soft
functions that obey the non-Abelian exponentiation theorem in this contribution. As a first
application of our method, we consider the $N$-jettiness observable and present numerical
results for the $1$-jettiness and $2$-jettiness hadron-collider soft functions to
next-to-next-to-leading order in the perturbative expansion.