This is now a reworked extended self-answer including all derivations.
First of all I would like to point out again, that the problem proposed in this OP is a natural generalization of finding closed forms for convergent sums to the field of divergent sums. Instead of asking for the value of a convergent sum and its possible closed form, we aks here for the asymptotic expansion which usually contains constants which tae the role of the convergent values.
The derivations will be explained later as soon as I find the time.
a) This first part can indeed by found, as suggested in a comment, by using the Euler-MacLaurin expansion for the partial sum in the form
$$\sum_{k=a}^b f(k) = \int_{a}^b f(x) \,dx + \frac{1}{2}(f(a) + f(b)) \\+ \sum_{j=1}^m \frac{B_{2j}}{(2j)!} \left(f^{(2j -1)}(b) - f^{(2j -1)}(a)\right) + R_m\tag{1}$$
where $f^{(k)}$(x) is the $k$-th derivative, $B_{j}$ is the $j$-th Bernoulli number ($B_2 = \frac{1}{6}$, $B_4 = -\frac{1}{30}$, ...), $m$ is a natural number, and $R_m$ is the remainder term to guarantee exact validity of the formula.
Setting $a=1$ and $b=n$, $f(x) = \sin(\sqrt{x})$, $f'(n)= \frac{\cos(\sqrt{n})}{2 \sqrt{n}}$, $f'''(n) = \frac{3 \sin \left(\sqrt{n}\right)}{8 n^2}-\frac{(n-3) \cos \left(\sqrt{n}\right)}{8 n^{5/2}}$ and $\int_0^n \sin(\sqrt{x})\,dx = 2 \sin \left(\sqrt{n}\right)-2 \sqrt{n} \cos \left(\sqrt{n}\right)$ we find from $(1)$
$$\sum_{k=0}^n \sin(\sqrt{k})=\\
\left\{2 \left(\sin \left(\sqrt{n}\right)-\sqrt{n} \cos \left(\sqrt{n}\right)\right)
+\frac{\sin \left(\sqrt{n}\right)}{2}
\\+\frac{\cos \left(\sqrt{n}\right)}{24 \sqrt{n}}
\\+\frac{1}{720} \left(\frac{\cos \left(\sqrt{n}\right)}{8 n^{3/2}}-\frac{3 \cos \left(\sqrt{n}\right)}{8 n^{5/2}}-\frac{3 \sin \left(\sqrt{n}\right)}{8 n^2}\right)+...\right\}\\
+\left[2 \cos (1)-2 \sin (1)+\frac{\sin (1)}{2}-\frac{\cos (1)}{24}+\frac{3 \sin (1)}{720\ 8}+\frac{\cos (1)}{720\ 4}+...\right]\tag{2} $$
Here we have separated two groups of terms: the curly brackets enclose all terms depending on $n$, the square brackets contain the constant terms.
Instead of explicitly calculating the sum of the constant terms we collect them all in a constant $c$. Our sum then becomes
$$s(n) = \sum_{k=0}^n \sin(\sqrt{k})= 2 \sin \left(\sqrt{n}\right)-2 \sqrt{n} \cos \left(\sqrt{n}\right)+\frac{1}{2}\sin(\sqrt{n}) \\+ \frac{1}{24}\left( \frac{\cos(\sqrt{n})}{ \sqrt{n}}\right)+c + O(n^{-\frac{3}{2}})\tag{3} $$
This is the formula which had to be proved. The final step consists in numerically determining the constant $c$ from inserting some appropriate value of $n$ into $(3)$ Here we have used $n=60$ to find a stable value of $c=-0.203569$.
b) I had just started a cumbersome method to find the constant $c$ but I dropped it immediately as I received the decisive hint in a comment of achille hui to use the Abel-Plana-formula (https://en.wikipedia.org/wiki/Abel%E2%80%93Plana_formula) which in its finite version(https://math.stackexchange.com/a/2541658/198592) reads.
$$\sum _{n=a}^b f(n)=\int_a^b f(x) \, dx+\frac{1}{2} (f(a)+f(b))\\
+i \int_{0}^{\infty } \frac{f(a+i t)-f(a-i t)-(f(b+i t)-f(b-i t))}{\exp (2 \pi t)-1} \, dt\tag{4}$$
Letting $a=1$, $b=n$ this gives for $f(x) = \sin(\sqrt{x})$
$$\sum _{n=1}^n f(n)=I_1(n)-I_1(1) + I_2(n)+c \tag{5}$$
where
$$I_1(n) = \frac{1}{2} \left(\sin \left(\sqrt{n}\right)\right)+2 \left(\sin \left(\sqrt{n}\right)-\sqrt{n} \cos \left(\sqrt{n}\right)\right)\tag{5a}$$
$$I_2(n) = 2 \int_0^{\infty } \frac{\cos \left(\frac{\sqrt{\sqrt{n^2+t^2}+n}}{\sqrt{2}}\right) \sinh \left(\frac{\sqrt{\sqrt{n^2+t^2}-n}}{\sqrt{2}}\right)}{\exp (2 \pi t)-1} \, dt\tag{5b}$$
and $c=I_2(0)$ is an integral representation of the constant defined in $a)$.
Let us first look at the integral $I_2(n)$.
For large $n$ the numerator of the integrand becomes
$$\sqrt{\frac{1}{n}} t \cos \left(\frac{\sqrt{\sqrt{n^2+t^2}+n}}{\sqrt{2}}\right)+O(n^{-\frac{3}{2}})$$
Noticing that due to the denominator contributions to the integral originate mainly from small $t$ we can assume $n >> t$ so that $\frac{\sqrt{\sqrt{n^2+t^2}+n}}{\sqrt{2}}$ can be replaced by $\sqrt{n}$ which gives
$$I_2(n) = \frac{\cos \left(\sqrt{n}\right)}{\sqrt{n}}
\int_0^{\infty } \frac{t}{\exp (2 \pi t)-1} \, dt
=\frac{\cos \left(\sqrt{n}\right)}{\sqrt{n}}+O(n^{-\frac{3}{2}}) \tag{6}$$
Hence we recover the expansion of a).
Next we turn to the constant, which is given by
$$c_i = -2 \int_0^\infty \frac{\sinh(\sqrt{\frac{t}{2}})\cos(\sqrt{\frac{t}{2}})}{e^{2 \pi t} - 1}\,dt\simeq -0.203569\tag{7}$$
The subscript $i$ indicates that we have an integral representation of $c$.
This integral most probably cannot be done in closed form. But we can derive interesting series representations.
Expanding the denominator in a power series as $\frac{1}{e^{2 \pi t} - 1} = \sum_{k=1}^\infty e^{-2 k \pi t} $ and doing the integrals results in the following formula
$$c_{s,1} = \frac{ \sqrt{2} }{4 \pi} \sum _{k=1}^{\infty } \frac{\cos \left(\frac{1}{8 \pi k}+\frac{\pi }{4}\right)}{ k^{3/2}} \simeq -0.203569\tag{8}$$
Expanding the $\cos$ in a power series and doing the $k$-sum results (after a lengthy manipulation involving binomial and the hypergeometric functions, see below) in this alternative formula for $c$
$$c_{s,2} = -\sum _{n=0}^{\infty } \frac{(-1)^{\frac{1}{4} \left(2 n-(-1)^n+1\right)} \zeta \left(\frac{1}{2} (2 n+3)\right)}{4 \pi (8 \pi )^n n!}\tag{9}\simeq -0.203569$$
Here $\zeta(s)$ is Riemann's zeta function. Notice that the sum in $(9)$ converges very quickly.
Derivation of $(9)$ from $(8)$
Expanding the $\cos(z) = \sum _{j=0}^{\infty } \frac{(-1)^j z^{2 j}}{(2 j)!}$ and expanding the binom a summand in $(8)$ becomes
$$s(j,k,m)=\frac{\sqrt{2} (-1)^j \left(\frac{\pi }{4}\right)^{2 j-m} \left(\frac{1}{8 \pi k}\right)^m \binom{2 j}{m}}{(4 \pi ) \left(k \sqrt{k}\right) (2 j)!}\tag{10}$$
Doing the $j$-sum (of the expandion of $\cos$) gives
$$s(k,m)= \sum_{j=1}^\infty s(j,k,m)= \frac{1}{k^{\frac{3}{2}}}2^{-m-\frac{3}{2}} \pi ^{-2 m-1} \binom{0}{m} \left(\frac{1}{k}\right)^m \, _1F_2\left(1;\frac{1}{2}-\frac{m}{2},1-\frac{m}{2};-\frac{\pi ^2}{64}\right)\tag{11}$$
where $_1F_2$ is a hypergeometric function.
The $k$-sum leads to
$$s(m)= \sum_{k=1}^\infty s(k,m)= 2^{-m-\frac{3}{2}} \pi ^{-2 m-1} \zeta \left(\frac{1}{2} (2 m+3)\right) b(m)\tag{12a}$$
where $\zeta$ is Riemann's zeta function and
$$b(m) = \binom{0}{m} \, _1F_2\left(1;\frac{1}{2}-\frac{m}{2},1-\frac{m}{2};-\frac{\pi ^2}{64}\right)\tag{12b}$$
Now the things are getting more complicated since, except for the case $m=0$, where $b(0) = \frac{\pi }{2 \sqrt{2}} J_{-\frac{1}{2}}\left(\frac{\pi }{4}\right)$, $b(m)$ is of the form $0*\infty$ for integer $m>0$.
In order to find limiting form of the product $b(m)$ for positive integer $m$ we write $\binom{0}{m} = \frac{\Gamma (1)}{\Gamma (m+1) \Gamma (1-m)}$ and use the series for the hypergeometric function
$$_1F_2\left(1;\frac{1}{2}-\frac{m}{2},1-\frac{m}{2};-\frac{\pi ^2}{64}\right) = \sum _{n=0}^{\infty } \frac{(1)_n z^n}{n! \left(\left(\frac{1-m}{2}\right)_n \left(1-\frac{m}{2}\right)_n\right)} = \sum _{n=0}^{\infty } \frac{z^n}{\left(\frac{1-m}{2}\right)_n \left(1-\frac{m}{2}\right)_n}\tag{13}$$
Here $(a)_n = \frac{\Gamma(a+n)}{\Gamma(n)}$ is the Pochhammer symbol and $z = -\frac{\pi ^2}{64}$.
Now we perform the limit under the $n$-sum and find for instance for $m\to 2$
$$\lim_{m\to 2} \, \frac{\Gamma (1) \Gamma \left(\frac{1-m}{2}\right) \Gamma \left(1-\frac{m}{2}\right) z^n}{(\Gamma (m+1) \Gamma (1-m)) \left(\Gamma \left(\frac{1-m}{2}+n\right) \Gamma \left(-\frac{m}{2}+n+1\right)\right)}= \frac{2^{2 n-1} z^n}{\Gamma (2 n-1)}$$
The first 6 Limits for $m=0..5$ are
$$\left\{\frac{4^n z^n}{\Gamma (2 n+1)},\frac{4^n z^n}{\Gamma (2 n)},\frac{2^{2 n-1} z^n}{\Gamma (2 n-1)},\frac{2^{2 n-1} z^n}{3 \Gamma (2 n-2)},\frac{2^{2 n-3} z^n}{3 \Gamma (2 n-3)},\frac{2^{2 n-3} z^n}{15 \Gamma (2 n-4)}\right\}\tag{14}$$
Doing the $n$-sum from $0$ to $\infty$ of the elements of this list and replacing $z$ we get
$$b(m= 0..5) = \left\{\frac{1}{\sqrt{2}},-\frac{\pi }{4 \sqrt{2}},-\frac{\pi ^2}{32 \sqrt{2}},\frac{\pi ^3}{384 \sqrt{2}},\frac{\pi ^4}{6144 \sqrt{2}},-\frac{\pi ^5}{122880 \sqrt{2}}\right\}\tag{15}$$
The denominators are found from https://oeis.org/A047053 to be $4^n n!$, and a formula for the sign pattern can be selected from https://oeis.org/A133872.
Inserting $(15)$ into $(12a)$ we have proved $(9)$.
