We have $$\left(\sqrt{x+1}-1\right)^4 = -4\sum_{m\geq 3}\left[\binom{1/2}{m-1}+2\binom{1/2}{m}\right]x^m $$ $$\left(\sqrt{1+\frac{1}{x}}-\sqrt{\frac{1}{x}}\right)^4 = -4\sum_{m\geq 1}\left[\binom{1/2}{m+1}+2\binom{1/2}{m+2}\right]x^{m} $$ hence: $$\begin{eqnarray*} \sum_{n\geq 1}\left(\sqrt{n+1}-\sqrt{n}\right)^4 &=& -4\sum_{m\geq 1}\left[\binom{1/2}{m+1}+2\binom{1/2}{m+2}\right]\,\zeta(m) \\&=&-4\int_{0}^{+\infty}\frac{dx}{e^x-1}\sum_{m\geq 1}\left[\binom{1/2}{m+1}+2\binom{1/2}{m+2}\right]\frac{x^{m-1}}{(m-1)!}\\&=&\int_{0}^{+\infty}\frac{2e^{-x/2}\left(x\,I_0\left(\frac{x}{2}\right)-4\,I_1\left(\frac{x}{2}\right)\right)}{x^2(e^x-1)}\,dx \\&=&2\int_{0}^{+\infty}\frac{I_2\left(x\right)\,dx}{x e^x (e^{x}-1)(e^x+1)}\end{eqnarray*}$$ where the last integrand function is pretty close to $\frac{1}{16}e^{-2x}$, from which it follows that the original series is pretty close to $\frac{1}{16}$. The following integral representation for the Bessel function $I_2$

$$ I_2(x) = \frac{x^2}{3\pi}\int_{0}^{\pi}\exp\left(x\cos\theta\right)\sin^4(\theta)\,d\theta $$ leads to: $$\begin{eqnarray*} \sum_{n\geq 1}\left(\sqrt{n+1}-\sqrt{n}\right)^4 &=& \frac{1}{6\pi}\int_{0}^{\pi}\psi'\left(\tfrac{3-\cos\theta}{2}\right)\sin^4(\theta)\,d\theta\\&=&-1+\frac{1}{6\pi}\int_{0}^{\pi}\psi'\left(\tfrac{1-\cos\theta}{2}\right)\sin^4(\theta)\,d\theta\\ &=&-1+\frac{\pi}{6}\int_{0}^{\pi/2}\frac{\sin^4(\theta)}{\sin^2\left(\pi\sin^2\frac{\theta}{2}\right)}\,d\theta\\ &=&\color{blue}{-1+\frac{\pi}{3}\int_{0}^{\pi/4}\frac{\sin^4(2\theta)\,d\theta}{\sin^2\left(\pi\sin^2\theta\right)}} \end{eqnarray*}$$ by the reflection formula for the trigamma function.
The blue integral can be written in the more symmetric form $$\boxed{\sum_{n\geq 0}\left(\sqrt{n+1}-\sqrt{n}\right)^4=\color{blue}{ \frac{4\pi}{3}\int_{0}^{1}\frac{x^{3/2}(1-x)^{3/2}}{\sin^2(\pi x)}\,dx}=\frac{\pi}{6}\int_{0}^{1}\frac{(1-x^2)^{3/2}\,dx}{\cos^2\frac{\pi x}{2}}.} $$


Let us tackle the case $s=3$ with a similar approach. We have $$(\sqrt{x+1}-1)^3 = \sum_{m\geq 3}\left[4\binom{1/2}{m}+\binom{1/2}{m-1}\right]x^m $$ $$\left(\sqrt{1+\frac{1}{x}}-\sqrt{\frac{1}{x}}\right)^3 = \sum_{m\geq 3}\left[4\binom{1/2}{m}+\binom{1/2}{m-1}\right]x^{m-3/2} $$ $$\begin{eqnarray*} \sum_{n\geq 1}\left(\sqrt{n+1}-\sqrt{n}\right)^3 &=& \sum_{m\geq 3}\left[4\binom{1/2}{m}+\binom{1/2}{m-1}\right]\zeta\left(m-\frac{3}{2}\right)\\&=&\sum_{m\geq 3}\left[4\binom{1/2}{m}+\binom{1/2}{m-1}\right]\frac{1}{\left(m-\frac{5}{2}\right)!}\int_{0}^{+\infty}\frac{x^{m-5/2}}{e^x-1}\,dx\\&=&\int_{0}^{+\infty}\frac{3 e^{-x} \left(2-2 e^x+x+e^x x\right)}{2 \sqrt{\pi } x^{5/2}(e^x-1)}\,dx\\&=&\frac{3}{\sqrt{\pi}}\int_{0}^{+\infty}\left[\frac{1}{2x^{3/2}e^x}-\frac{1}{x^{5/2}e^x}+\frac{1}{x^{5/2}e^x(e^x-1)}\right]\,dx\end{eqnarray*}$$ hence $$\boxed{\sum_{n\geq 0}\left(\sqrt{n+1}-\sqrt{n}\right)^3 = \color{blue}{\tfrac{3}{2\pi}\,\zeta\left(\tfrac{3}{2}\right)}}$$ just follows from integration by parts, Frullani's Theorem and the integral representation for the $\zeta$ function. The same approach allows an explicit evaluation in terms of the $\zeta$ function for any odd value of $s$. For instance:

$$\boxed{ \sum_{n\geq 0}\left(\sqrt{n+1}-\sqrt{n}\right)^5 = \color{blue}{\tfrac{15}{2\pi^2}\,\zeta\left(\tfrac{5}{2}\right)},\qquad \sum_{n\geq 0}\left(\sqrt{n+1}-\sqrt{n}\right)^7 = \color{blue}{\tfrac{7}{2\pi}\,\zeta\left(\tfrac{3}{2}\right)-\tfrac{105}{2\pi^3}\,\zeta\left(\tfrac{7}{2}\right)}\\\qquad\qquad\qquad\sum_{n\geq 0}\left(\sqrt{n+1}-\sqrt{n}\right)^9=\color{blue}{\tfrac{90}{2\pi^2}\zeta\left(\tfrac{5}{2}\right)-\tfrac{945}{2\pi^4}\zeta\left(\tfrac{9}{2}\right)}.}$$


2018 Update. This technique turns out to be pretty flexible. By following the same steps, i.e.

  1. Expanding $g(n)$ as a power series in $\frac{1}{n}$;
  2. Exploiting the integral representation for $\zeta(s)$ and switching $\sum,\int$;
  3. Exploiting the integral representation for Bessel functions of the first kind and switching $\int_{0}^{+\infty},\int_{0}^{\pi}$;
  4. Employing integration by parts

one may also prove $$ \sum_{n\geq 1}\left(\sqrt{n^2+1}-n\right)^2 = 2\int_{0}^{+\infty}\frac{J_2(x)\,dx}{x(e^x-1)}=\frac{2}{\pi}\int_{0}^{1}\sqrt{1-x^2}\left(\pi x\coth(\pi x)-1\right)\,dx $$ and similar identities. This is pretty reminiscent of the Voronoi summation method.


It's straightforward to note the case $\phi(3)$ (and similarly for the case $\phi(5)$) is

$$1+\lim_{s\to-3/2}\left(4\sum_{n=1}^{\infty} \frac{1}{(n+1)^s}-3\sum_{n=1}^{\infty} \frac{1}{(n+1)^{s+1}}-4\sum_{n=1}^{\infty} \frac{1}{n^s}-3\sum_{n=1}^{\infty} \frac{1}{n^{s+1}}\right)$$ $$=-6\lim_{s\to-3/2}\zeta(1+s)=-6\zeta(-1/2)=3/(2\pi) \zeta(3/2),$$ where the last equality follows directly from the Riemann functional equation here.

Q.E.D.