I tried to answer by proposing a simple model. The model didn't work as well as I hoped, but the general procedure might serve as an answer.
A harmonic model
As a first approximation, you might consider a simple model for the energy of a single water molecule as a function of its angle. One possibility might be a harmonic function
$$E(\theta) = \frac{k \theta^2}{2}$$
where $k$ is an unknown parameter and $\theta$ is the angle. This is very simple and leaves out many complexities, of course.
If you had $E(\theta)$, the energy of a single water molecule at a particular angle, and if you assume an independent particle model (which may break down for a liquid), you could apply Boltzmann weighting
$$p(\theta) = \frac{1}{U}\exp\left(-\frac{E(\theta)}{R T}\right)$$
where $p(\theta)$ is the probability of finding a water molecule at a particular angle and $U = \int_0^\pi d\theta \exp\left(-\frac{E(\theta)}{R T}\right)$ is a normalization factor.
If you had a probability distribution $p(\theta)$ describing the spread of angles at a certain temperature, getting the mean and standard deviation is easy:
$$\langle\theta\rangle = \int_0^\pi d\theta\ \theta\ p(\theta)\\
\langle\theta^2\rangle = \int_0^\pi d\theta\ \theta^2\ p(\theta)\\
\sigma_\theta^2 = \int_0^\pi d\theta\ \left(\theta - \langle\theta\rangle\right)^2\ p(\theta) = \langle\theta^2\rangle - \langle\theta\rangle^2$$
So this would be how to get $\sigma_\theta^2$.
Below is my numerical attempt.
Beware: this model is too simple.
Numerical results
I did some calculations (mostly with SymPy) with the above and I came up with the following:
$$p(\theta) = \frac{1}{U} \exp\left(-\frac{k \theta^2}{2 R T}\right)\\
U = \sqrt{\frac{\pi R T}{2 k}} erf \left(\frac{\pi}{2} \sqrt{\frac{2 k}{R T}} \right)$$
$$\langle\theta\rangle = \frac{k}{R T U}\left(\exp\left( \frac{\pi^2 k}{2 R T} \right) -
1 \right) \exp\left( -\frac{\pi^2 k}{2 R T} \right)\\
\sigma_\theta^2 = \frac{R T}{k} - \left[ \frac{\pi}{U} \exp\left( -\frac{\pi^2 k}{2 R T} \right) +
\left(\frac{RT}{kU}\right)^2 \left(\exp\left( \frac{\pi^2 k}{2 R T} \right) -
1 \right)^2 \exp\left( -\frac{\pi^2 k}{R T} \right)
\right]$$
By substituting constants and temperature ($\pi = 3.14159265$, $R = 8.314459848$ J K$^{-1}$mol$^{-1}$ and $T = 293.15$ K) and setting the average above to as close as possible to 104.45° and solving numerically the complete least squares problem, I found an approximate solution for $k$: $384.05105591$ J/rad².
Sadly, this value is associated with an average angle of around 96°, which means the model is very poor.
Using this best $k$ gives us a $\sqrt{\sigma_\theta^2}$ of
$$\sqrt{0.508770928037803} = 0.713281801280393 \text{ rad}$$
which is around 41°. Again, a too simplistic model: the best prediction gives angles as 96°$\pm$41°!
A bad model
By expanding the average in series, I found the following
$$\langle\theta\rangle \approx 1.570796325 + 0.000530047067361739 k + O(k^2)$$
The second order term above was $- 7.86976333117302\times 10^{-7} k^2$, so I dropped it.
Observe the first term is pretty much 90° in radians. All other terms decay too quickly.
By doing the same expansion for the standard deviation the following comes out:
$$\sigma_\theta^2 \approx - \frac{4.54747350886464 \times 10^{-13}}{k} + 0.411233515772252 + 0.000388544793228128 k + O(k^2)$$
The next term in the series above was $- 4.34803406790969 \times 10^{-7} k^2$ so I dropped it. Again, a (wrong) zeroth order term dominates.
