I found a partial answer inspired by the paper "Approximate fixed point iteration with an application to infinite horizon Markov decision processes." The gist is that
$$\limsup_n E[|S_n - r^2|] \leq C \sigma,$$
for some constant $C$, where $r = E[R_n]$ and $\sigma^2 = E[(R_n - r)^2]$. So $S_n$ may not actually converge to $r^2$, but it is on average not that far.
Say that your rewards are centered around an average value $r > 1$. That is, $R_n = r + \epsilon_n$, where $\epsilon_n$ is a sequence of identically distributed random variables with mean $0$ and variance $\sigma^2$. Then our equation reads $$S_{n + 1} = S_n + r + \epsilon_{n + 1} - \sqrt{S_n}.$$
If we let $T(x) = x + r - \sqrt{x}$, then we can write $$S_{n + 1} = T(S_n) + \epsilon_{n + 1}.$$ This shows that we can think of $S_n$ as the result of iterating the map $T(x)$ with some "noise." Since $T(r^2) = r^2$, an application of the mean value theorem gives
$$|T(x) - r^2| \leq |x - r^2| \left(1 - \frac{1}{2 \sqrt{z}} \right)$$
for some $z$ between $x$ and $r^2$. Thus $T$ is a contraction map on any suitably bounded interval. Let's write $|T(x) - r^2| \leq c |x - r^2|$ for some $0 < c < 1$.
All of this together gives
\begin{align*}
|S_{n + 1} - r^2| &\leq |T(S_n) - r^2| + |\epsilon_{n + 1}| \\
&\leq c |S_n - r^2| + |\epsilon_{n + 1}|.
\end{align*}
If we set $E_n = E[|S_n - r^2|]$, then taking expectations of this above equation yields
$$E_{n + 1} \leq c E_n + E[|\epsilon_{n + 1}|] \leq c E_n + \sigma.$$
Repeating this argument gives
$$E_n \leq c^n E_0 + \sigma \sum_{k = 0}^{n - 1} c^{n-k} = c^n E_0 + \frac{c \sigma}{1 - c} (1 - c^n).$$
The right-hand side goes to $c \sigma / (1 - c) = C \sigma$ as $n \to \infty$.
To be really rigorous we need to prove that $S_n$ is almost surely contained in an interval where $T$ is a contraction. I'm not sure how to do that, but this seems like a reasonable first sketch of what's going on.