If $\;z\in(1,\infty),\;$ then the function
$$f_n(z)=\dfrac\pi4\,(z^2-1)\,z^{\large-\frac{n+4}2}\cot{\pi\dfrac{(z-1)^2}{4z}}\tag1$$
has the poles in the points
$$z_k = 2k+1 +2\sqrt{k^2+k} =\left(\sqrt{k+1}+\sqrt k\right)^2,\quad k\in\mathbb N,\tag2$$
wherein
$$\underset{\large {z_k}\small^\mathstrut}{\operatorname{Res}}f_n(z)
=\dfrac\pi4\,(z^2_k-1)z_k^{\large-\frac{n+4}2}\,{\cos\pi\dfrac{(z-1)^2}{4z}}
\lim\limits_{\large z\to {z_k}\Large\mathstrut}\,\dfrac{z-z_k}{\sin\pi\dfrac{(z-1)^2}{4z}}
$$$$
=\dfrac\pi4\,(z^2_k-1)z_k^{\large-\frac{n+4}2}\,{\cos\pi\dfrac{(z-1)^2}{4z}}
\lim\limits_{\large z\to {z_k}\Large\mathstrut}\,
\dfrac1{\dfrac\pi4\left(1-\dfrac1{z^2}\right)\cos\pi\dfrac{(z-1)^2}{4z}}
= z_k^{\large-\frac n2},$$
$$\underset{\large {z_k}\small^\mathstrut}{\operatorname{Res}}f_n(z)
= \left(\sqrt{k+1} - \sqrt k\right)^n.\tag3$$
Choosen the arbitrary infinity contour $\;C\;$ with the points $\;z_1, z_2,\dots,\;$
one can get
$$S_n=\sum\limits_{k=0}^\infty\left(\sqrt{k+1} - \sqrt k\right)^n
= 1+\dfrac1{2\pi i}\oint\limits_C f_n(z)\text dz\tag4$$
In particular, for the contours
$$C =\left\{\dfrac{1}{1-t^2}+p+32it,\;t\in(-1,1),\;p\in[1,4]\right\}$$

easily to get the integral representation in the form of
$$S_n= 1+\dfrac1{2\pi i}\int\limits_{-1}^1 f_n\left(\dfrac{1}{1-t^2}+p +32it\right)\left(\dfrac{2t}{(1-t^2)^2}+32i\right)\text dt.\tag5$$
On the other hand,
$$\left(\sqrt{k+1}-\sqrt k\right)^n = \left(2k+1+\sqrt{(2k+1)^2-1}\right)^{\large-\frac n2} $$
$$= (2k+1)^{\large-\frac n2}\left(1+\sqrt{1-\dfrac1{(2k+1)^2}}\right)^{\large-\frac n2}$$
$$= (4k+2)^{\large-\frac n2}\left(1-\dfrac1{(4k+2)^2}-\dfrac1{(4k+2)^4}-\dots\right)^{\large-\frac n2}$$
$$= (4k+2)^{\large-\frac n2}\left(1+\dfrac n{2(4k+2)^2}+\dfrac n{2(4k+2)^4}
+\dfrac{n(n+2)}{8(4k+2)^4}+\dots\right),$$
$$\begin{align}
&S_n\approx \sum\limits_{k=0}^{s-1}\left(\sqrt{k+1} + \sqrt k\right)^{-n} + 2^{-n}\\[4pt]
&\times\left(\zeta\left(\dfrac n2, s+\dfrac12\right)
+\dfrac n{32}\zeta\left(\dfrac n2+2, s+\dfrac12\right)
+\dfrac{n(n+6)}{2048}\zeta\left(\dfrac n2+4, s+\dfrac12\right) \right), \end{align}\tag6$$
where $\;\zeta(z,a)\;$ is the generalized Riemann zeta function.
Numeric calculations of the given sum via $(5)$ (with the same results for all values of $\;p\;$ from domain) and via $(6)$ for $\;s=35\;$ by Wolfram Alpha are shown in the table $(7).$
\begin{vmatrix}
n & (5) & (6)\\
3 & 1.24731\,73498\,64127\,39610\,384 & 1.24731\,73498\,64124\,19\\
4 & 1.06038\,21964\,75434\,10179\,32212 & 1.06038\,21964\,75433\,7130\\
5 & 1.01940\,80755\,93321\,43798\,17195 & 1.01940\,80755\,93321\,3941\\
6 & 1.00693\,53143\,70406\,29845\,53855\,17 & 1.00693\,53143\,70406\,29373\,3\\
7 & 1.00261\,51034\,44489\,88275\,39751\,41 & 1.00261\,51034\,44489\,88226\,2\\
8 & 1.00101\,71884\,92355\,62561\,09186\,051 & 1.00101\,71884\,92355\,62556\,09\\
9 & 1.00040\,33950\,87623\,01165\,83908\,220 & 1.00040\,33950\,87623\,01165\,34\\
10 & 1.00016\,20240\,77476\,98571\,45268\,798 & 1.00016\,20240\,77476\,98571\,404\\
11 & 1.00006\,56411\,98162\,21005\,94527\,8947 & 1.00006\,56411\,98162\,21005\,9406\\
12 & 1.00002\,67536\,86676\,80665\,73794\,1985 & 1.00002\,67536\,86676\,80665\,73732\tag7
\end{vmatrix}
Therefore, the calculations confirm obtained formulas $(4)-(5).$