Is your x,y,z data set gridded? If not, grid it, perhaps with gridfit. With your gridded data set (I'll distinguish with capital letters), you'll have some regular grid of Z values corresponding to X and Y values. If X and Y are at, say 100 meter spacing, the area of each grid box is 100x100 meters = 10,000 square meters. A matrix of volume of water in each grid box is then given by
VolumeMap = (SurfaceElevation - Z)*(100*100).
And total volume of water in cubic meters is given by
sum(VolumeMap(:))
