Learn smoothly varying mean and variance of a variable over a 2d domain

For a problem which I am working on at the moment, I'm interested in learning how the mean and variance of some response variable y changes with two independent variables x1 and x2 - i.e. for each coordinate in (x1, x2)-space I wish to have an estimate for $\mu_y$ and $\sigma_y$ in order to be able to approximately standardise new observations as they arrive.

I have enough domain knowledge to expect both the mean and variance of y to vary smoothly across this space (plus I can likely place bounds on the gradient at any point), which I'd like to use to be able to deal with previously-unseen x1, x2 pairs as they arise.

What techniques would you suggest for such a problem?