I will make the first attempt, though my resulting bound can certainly be improved.

I will do the $N\times N$ Stirling matrix of the second kind. We have $S_2(n,k) = \left\{n\atop k\right\}$ for $0\leq n,k\leq N-1$, which satisfies the following recurrence relation:

$$\left\{n+1\atop k\right\} = k\left\{n\atop k\right\} + \left\{n\atop k-1\right\},\qquad n\geq 0, \quad k>0.$$

Therefore, the matrix $S_2$ satisfies the Sylvester matrix equation:

$$\underbrace{\begin{bmatrix}0 & 1 \\ & & \ddots \\ & & & 1 \\ -1 \end{bmatrix}}_{=A} S_2 - S_2\underbrace{\begin{bmatrix} 0 & 1 \\ & 1 & \ddots \\ & & \ddots & 1 \\ & & & N-1 \end{bmatrix}}_{=B} = \begin{bmatrix}0 & \ldots & 0 & 0 \\\vdots & \ddots & \vdots & \vdots \\0 & \ldots & 0 & 0 \\ \times & \ldots & \times & \times \end{bmatrix},$$
where the $\times$'s denote nonzero entries. We say that $S_2$ has an $(A,B)$-displacement rank of 1.

From now on assume that $N$ is an even integer. The same idea, with different details, works when $N$ is an odd integer.

Since the rhs of the above equation is of rank 1, the eigenvalues of $A$ lie at (shifted) roots-of-unity, and the eigenvalues of $B$ are in $\{0,\ldots,N-1\}$ we find that, see Corollary 2.2 of paper
$$\sigma_{2k+1}(S_2) \leq \|V\|_2\|V^{-1}\|_2 Z_{2k}(\{0,\ldots,N-1\},F_+\cup F_-)\|S_2\|_2,$$
where $V$ is the eigenvector matrix for $B$, i.e., $B = V\Lambda_B V^{-1}$ and $\|\cdot\|_2$ is the spectral norm. Here, $Z_{k}(E,F)$ denotes a Zolotarev number (see Section 2 of paper) and

$$F_+ = \left\{e^{it}: t\in [\tfrac{\pi}{N},\pi-\tfrac{\pi}{N}\right\}, \quad F_- = \left\{e^{it}: t\in [-\pi + \tfrac{\pi}{N},-\tfrac{\pi}{N}\right\}.$$

I do not have bounds on $Z_{2k}(\{0,\ldots,N-1\},F_+\cup F_-)$ at hand and would have to work them out. Instead, to use a previous result, I will use the fact that

$$Z_{2k}(\{0,\ldots,N-1\},F_+\cup F_-)\leq Z_{2k}(\mathbb{R},F_+\cup F_-)$$

since $\{0,\ldots,N-1\}\subset \mathbb{R}$.

Therefore, when $N$ is an even integer, we obtain (see Lemma 5.1 of paper):

$$\sigma_{2k+1}(S_2) \leq 4\|V\|_2\|V^{-1}\|_2\left[\exp\left(\frac{\pi^2}{4\log(4N/\pi)}\right)\right]^{-k}\|S_2\|_2.$$

This explains the rapid decay of the singular values of $S_2$, provided $\|V\|_2\|V^{-1}\|_2$ does not grow rapidly with $N$. With a little more work, I believe that one can show that $\|V\|_2\|V^{-1}\|_2<8$ for any $N$ (this appears to be true numerically for $1<N<1000$). If so, an explicit bound on the singular values is (for $N$ even):

$$\sigma_{2k+1}(S_2) \leq 32\left[\exp\left(\frac{\pi^2}{4\log(4N/\pi)}\right)\right]^{-k}\|S_2\|_2, \qquad 0\leq 2k\leq N-1.$$

DISCLAIMER: This bound explains the geometric decay of the singular values of $S_2$ but the geometric rate is not tight. One would have to directly bound $Z_{2k}(\{0,\ldots,N-1\},F_+\cup F_-)$ to get a tighter bound. I do not know how to explain the two different rates of decay of the singular values yet.