I'm going to solve this by calculation. I'm sure there's some clever statistics way to figure this out -- if you know, please make a post and I'll probably accept it :p. Until then, this will have to do.
Problem statement
We know $\mu_1=\mu(S_1)$, $\sigma_1=\sigma(S_1)$, and $N_1=\lvert S_1\rvert$, as well as the equivalent statistics on $S_3$: $\mu_3$, $\sigma_3$, and $N_3$. We also know $N_2$ and wish to find estimates for $\mu_2$ and $\sigma_2$ along with error bounds.
Statistical identities on disjoint unions
For $X\cap Y=\emptyset$ we'll make frequent use of the identities
$$
\mu_{X\cup Y}=\frac{N_X\mu_X+N_Y\mu_Y}{N_X+N_Y}
$$
and
$$
(N_X+N_Y)\sigma_{X\cup Y}^2=N_X\sigma_X^2+N_Y\sigma_Y^2+\frac{N_XN_Y}{N_X+N_Y}(\mu_X-\mu_Y)^2.
$$
S3 minus S1
Using these formula, we can immediately find exactI mean and variance for $S_3\setminus S_1$, namely
$$
\mu_{3\setminus1}=\frac{N_3\mu_3-N_1\mu_1}{N_3-N_1}
$$
and
$$
(N_3-N_1)\sigma_{3\setminus1}^2=N_3\sigma_3^2-N_1\sigma_1^2-\frac{N_1N_3}{N_3-N_1}(\mu_3-\mu_1)^2.
$$
S2 minus S1
If we can use this information to find good estimates for $\mu_{2\setminus1}$ and $\sigma_{2\setminus1}$, the mean and standard deviation of $S_2\setminus S_1$, then we can use the union formulae to make estimates for $\mu_2$ and $\sigma_2$. Intuitively we expect that $\mu_{2\setminus1}\approx\mu_{3\setminus1}$ and $\sigma_{2\setminus1}\approx\sigma_{3\setminus1}$ are reasonable estimates on that subset, but we would need to justify this, at least finding upper and lower bounds to determine how good (or bad) these estimates are.
Now we have by identity an equation we'll repeatedly refer to as (A)
$$
N_{3\setminus1}\sigma_{3\setminus1}^2=N_{3\setminus2}\sigma_{3\setminus2}^2+N_{2\setminus1}\sigma_{2\setminus1}^2+\frac{N_{3\setminus2}N_{2\setminus1}}{N_{3\setminus1}}(\mu_{3\setminus2}-\mu_{2\setminus1})^2
$$
where we've abbreviated $N_{a\setminus b}:=\lvert{S_a\setminus S_b\rvert}=N_a-N_b$. If we wish to maximize $\sigma_{2\setminus1}$ it suffices to note that we could choose a distribution of $S_3\setminus S_1$ so that all the varianceII is in $S_2\setminus S_1$ and the means $\mu_{3\setminus2}-\mu_{2\setminus1}=0$. In this case, all the terms of (A) become zero except the LHS and the $\sigma_{2\setminus1}$ term, yielding
$$
N_{3\setminus1}\sigma_{3\setminus1}^2=N_{2\setminus1}\sigma_{2\setminus1}^2,
$$
implying that $\sigma_{2\setminus1}^2\leq\frac{N_{3\setminus1}}{N_{2\setminus1}}\sigma_{3\setminus1}^2$. By the same reasoning we could choose a different distribution so that all the variance is in $S_3\setminus S_2$, in which case $\sigma_{2\setminus1}^2=0$, obviously a lower bound. So we have these bounds on $\sigma_{2\setminus1}$:
$$
0\leq\sigma_{2\setminus1}\leq\sqrt{\frac{N_3-N_1}{N_2-N_1}}\sigma_{3\setminus1}.
$$
To find the maximum difference of the mean $\mu_{2\setminus1}$ from our proposed estimate $\mu_{2\setminus1}$ a good first step is to note we can maximize $(\mu_{3\setminus2}-\mu_{2\setminus1})^2$ in (A) by choosing a distribution where $\sigma_{3\setminus2}=\sigma_{2\setminus1}=0$. In this case all the variance in $S_3\setminus S_1$ is due to the wide disparity between the means on $S_3\setminus S_2$ and $S_2\setminus S_1$. So from (A) we get
$$
\lvert\mu_{3\setminus2}-\mu_{2\setminus1}\rvert\leq\frac{N_{3\setminus1}}{\sqrt{N_{3\setminus2}N_{2\setminus1}}}\sigma_{3\setminus1}.
$$
If we rearrange this along with our mean union identity at the top, we get two equations with two unknowns:
$$
\begin{align}
N_{3\setminus2}(\mu_{3\setminus2}-\mu_{3\setminus1}) +N_{2\setminus1}(\mu_{2\setminus1}-\mu_{3\setminus1}) &=0 \\
(\mu_{3\setminus2}-\mu_{3\setminus1}) -(\mu_{2\setminus1}-\mu_{3\setminus1}) &={\frac{N_{3\setminus1}}{\sqrt{N_{3\setminus2}N_{2\setminus1}}}}\sigma_{3\setminus1}
\end{align}
$$
Solving this gives us our bounded error:
$$
\left\lvert\mu_{2\setminus1}-\mu_{3\setminus1}\right\rvert\leq\sqrt{\frac{N_3-N_2}{N_2-N_1}}\sigma_{3\setminus1}.
$$
So now all the pieces are in place, and we need to simply combine what we got on $S_2\setminus S_1$ with what we already know for $S_1$ and we will have our estimate and error bounds.
S2
Now the estimated mean and its bounds on $S_2$ can readily be found:
$$
\begin{align}
\mu_2&=\frac{N_1\mu_1+(N_2-N_1)\left(\mu_{3\setminus1}\pm\sqrt{\frac{N_3-N_2}{N_2-N_1}}\sigma_{3\setminus1}\right)}{N_2}\\
&=\frac{N_1\mu_1+(N_2-N_1)\mu_{3\setminus1}}{N_2}\pm\frac{\sqrt{(N_2-N_1)(N_3-N_2)}}{N_2}\sigma_{3\setminus1}\\
\end{align}
$$
These bounds are exact and this estimate is quite good when $N_3-N_1\ll N_2$ and $\sigma_{3\setminus1}$ is small.
Solving for $\sigma_2$, we find that
$$
\sigma_2^2=\frac{N_{2\setminus1}}{N_2}\sigma_{2\setminus1}^2+\frac{N_1}{N_2}\sigma_1^2+\frac{N_{2\setminus1}N_1}{N_2}(\mu_{2\setminus1}-\mu_1)^2
$$
We won't find the exact error here because $\mu_{2\setminus1}$ and $\sigma_{2\setminus1}$ aren't independent, but the $v$ term gives rise to an error of order $\frac{N_{3\setminus1}}{N_2}\sigma_{3\setminus1}^2$ and the squared difference of means gives rise to the (small if the mean error is small) error term
$$
\frac{N_{2\setminus1}^2N_{3\setminus2}N_1}{N_2^4}\sigma_{3\setminus1}
$$
and the potentially more concerning error term
$$2\left\lvert\mu_{3\setminus1}-\mu_1\right\rvert\frac{N_{2\setminus1}N_1\sqrt{N_{2\setminus1}N_{3\setminus2}}}{N_2^3}\sigma_{3\setminus1}.
$$
ITheoretically we can find this exactly. In practice there are precision issues with this calculation. In my particular case I actually have these values already, but if you don't, beware!
II There's a actually special cases here where $N_{3\setminus2}=1$ or $N_{3\setminus2}=1$ since one element sets don't have any variance. But the bounds derived still hold -- we could just have made them tighter for these cases.
