This probability can be calculated exactly, and indeed it approaches $1/2$ rather quickly — more precisely, for each $p$ it approaches the fraction $(p-1)/(2p)$ of quadratic residues $\mod p$$\bmod p$. This can be proved by elementary means, but perhaps the nicest way to think about it is that if you choose $n$ numbers $a_i$ independently and sum $a_i^2 \bmod p$, the resulting distribution is the $n$-th convolution power of the distribution of a random single square — so its discrete Fourier transform is the $n$-th power of the D.F.T., call it $\gamma$, of the distribution of $a^2 \bmod p$. For this purpose $\gamma$ is normalized so $\gamma(0)=1$. Then for $\gamma(k) = (k/p) \gamma(1)$$k \neq 0$ we have (where$\gamma(k) = (k/p) \gamma(1)$ [where $(\cdot/p)$ is the Legendre symbol)symbol], and $p \gamma(1) = \sum_{a\bmod p} e^{2\pi i a^2/p}$, which $$ p \gamma(1) = \sum_{a \bmod p} \exp(2\pi i a^2/p), $$ which is a Gauss sum and is thus a square root of $\pm p$. It follows that $|\gamma(k)| = p^{-1/2}$, from which we soon see that each value of the convolution approaches $1/p$ at the exponential rate $p^{-n/2}$, and the probability you asked for approaches $(p-1)/(2p)$ at the exponentialsame rate $p^{-n/2}$.
As noted above, this result, and indeed the exact probability, can be obtained by elementary termsmeans, yielding a (known but not well-known) alternative proof of Quadratic Reciprocity(!). But that's probably too far afield for the present purpose.