This sounds like a hierarchical model (or multi-level model), but not explained in those terms. Firstly, note that a Poisson "residual" is not really the same thing as you would find in the case of a normal error term. The reason is that a poisson random variable cannot be separated out easily into a "systematic part" and "random part". One standard way is to use the "deviance residuals", defined by:
$$d_{i}=sgn(Y_{i}-\hat{Y}_{i})\sqrt{2Y_{i}\log\left(\frac{Y_{i}}{\hat{Y}_{i}}\right)-2(Y_{i}-\hat{Y}_{i})}$$
These quantities are analogues to the residuals in OLS regression, and can be easily obtained from residuals(glm.object, "deviance") function in R. You can show that the sum of squared deviance residuals has an approximate chi-square distribution, but I wouldn't know how to figure out the exact distribution.
An alternative is to add further structure to the Poisson GLM by introducing "random effects" into the model equation. For the Poisson model with canoncial log-link function we have $\log(\mu_i)=x_{i}^{T}\beta$. In your case $x_{i}^{T}$ is a row of a "design matrix" (all zeros and ones) which picks out the relevant beta coefficient for your site effect and year effect. This addition of random effects basically says that the above relationship only holds "on the average" rather than deterministically (can think of this as a "softer" link function). So that we actually have $E\left[\log(\mu_i)\right]=x_{i}^{T}\beta$. The "actual link" is given by:
$$\log(\mu_i)=x_{i}^{T}\beta+z_{i}^{T}u_{i}$$
Where $z_{i}^{T}$ is a known set of covariates, and $u_i$ is a random effect for the $i$th observation with $E[u_{i}]=0$. Now if you were to set $z_{i}^{T}=1$ then you have a good definition of the residuals, given by the random effects $u_i$. All that is needed to do now is to assign a probability distribution to the random effects and do the maths/coding work to either give point estimates or posterior probability distributions for $u_{i}$ to use in your subsequent model. These are often called Generalised linear mixed models, and $u_i$ is often chosen to have a normal distribution.
The nlme package should be able to do this kind of modelling
I've been estimating continuous positive outcome Poisson regressions with the Huber/White/Sandwich linearized estimator of variance fairly frequently. However, that's not a particularly good reason to do anything, so here are some actual references.
From the theory side, $y$ does not need to be an integer for for the estimator based on the Poisson likelihood function to be consistent. This is shown in Gourieroux, Monfort and Trognon (1984). This is called Poisson PMLE or QMLE, for Pseudo/Quasi Maximum Likelihood.
There's also some encouraging simulation evidence from Santos Silva and Tenreyro (2006), where the Poisson comes in best-in-show. It also does well in a simulation with lots of zeros in the outcome. You can also easily do your own simulation to convince yourself that this works in your snowflake case.
Finally, you can also use a GLM with a log link function and Poisson family. This yields identical results and placates the count-data-only knee jerk reactions.
References Without Ungated Links:
Gourieroux, C., A. Monfort and A. Trognon (1984). “Pseudo Maximum Likelihood Methods: Applications to Poisson Models,” Econometrica, 52, 701-720.
Best Answer
The key assumption of a generalized linear model that's relevant here is the relationship between the variance and mean of the response, given the values of the predictors. When you specify a Poisson distribution, what this implies is that you are assuming the conditional variance is equal to the conditional mean.* The actual shape of the distribution doesn't matter as much: it could be Poisson, or gamma, or normal, or anything else as long as that mean-variance relationship holds.
* You can relax the assumption that the variance equals the mean to one of proportionality, and still usually get good results.