Performing a Spatial Join will do this. Right click on the point layer and choose "Joins and Relates > "Join". In the Join Data dialog box, choose "Join data from another layer based on spatial location" in the drop down. Then choose the polygon layer you want joined. Then choose the radio button that says "is closest to it". (The selections are a little different if you're joining points to lines)
This function is also available in ArcToolbox, which is supposed to give better performance with large datasets, and provides some extra functionality.
This method is available in all versions of ArcGIS, and at all license levels.
Perform a hydrologic analysis on your data. Taking your first step of water bodies as a raster, you can then use that as a sinks raster. I'll specify the rest of this analysis in terms of GRASS as you mentioned that's the GIS system you're using:
Set up r.watershed
(documentation) with the elevation data layer you'd like to analyze and the sinks raster you generated in the first step:
r.watershed elev=input_dem depression=input_sinks basin=output_basins \
stream=output_streams threshold=1000
Where threshold is something appropriate for the scale of your data-- this should give you a map of each watershed: if you subtract that elevation from all the cells in that basin, you should get the vertical distance to the nearest water. You may need to iterate over smaller regions than your full raster to get good performance. You might find these tutorials (1, 2) also helpful in understanding how to use the command.
Mike helpfully mentioned an addon called r.watershed.distance which can be used to calculate this all in one go:
r.stream.distance -o dir=dirs stream=streams dem=elev \
distance=distance_outlets elevation=elevation_outlets
Which will result in an output close to what you're interested in:
This example was taken from the R.stream.* page, you can download the extension itself from the addon page.
Best Answer
If you are working with GRASS, then you can use the module v.distance. When you import the 2kmX2km grid into grass, it automatically creates centroids for each grid cell (that's how the vector design in GRASS works). So you would first add a column in the grids table to hold the distance value, then run v.distance to update that column:
(Assuming two GRASS vector maps, 'water' and 'grid')
View the results: