The function $A\rightarrow S=\sqrt{A}$ is defined and differentiable on the set of SPD matrices. Let $K=DS_A(H)$ be the derivative of $S$ in $A$, where $H$ is a variable SYMMETRIC matrix. Here $SS=A$ implies $KS+SK=H$, a Sylvester equation in the unknown $K$. We may assume $A=diag((\lambda_i)_i)$ where $\lambda_i>0$. Then $S=diag((\sqrt{\lambda_i})_i)$. If $H=[h_{i,j}],K=[k_{i,j}]$, then, by an easy identification, we obtain $k_{i,j}=\dfrac{h_{i,j}}{\sqrt{\lambda_i}+\sqrt{\lambda_j}}$. Of course, if $n=1$, we find the usual derivative $h_{1,1}\rightarrow \dfrac{h_{1,1}}{2\sqrt{\lambda_1}}$.
EDIT 1. About the half-vectorization operator, we can store half of matrix $K$ because it is symmetric as the matrix $H$.
EDIT 2. Another form of $K$ is $\int_0^{\infty}e^{-tS}He^{-tS}dt$. That implies that if $H$ is a small symmetric matrix, then $\sqrt{A+H}\approx \sqrt{A}+\int_0^{\infty}e^{-t\sqrt{A}}He^{-t\sqrt{A}}dt$.
EDIT 3. Proof of the above result. The integral converges (easy) and it suffices to prove that $KS+SK=H$ (the solution in $K$ of this equation is unique). One has $KS+SK=\int_0^{+\infty}e^{-tS}He^{-tS}S+Se^{-tS}He^{-tS}dt=-\int_0^{+\infty}(e^{-tS}He^{-tS})'dt=H.$