I think this is best handled in the interpolation process rather than on the back end of the interpolation process.
My suggestion would be to use 3D Analyst to create a terrain dataset from your bare earth LiDAR points.
- Create the terrain (3D Analyst>Terrain Management>Create Terrain)
- Add your LiDAR points (3D Analyst>Terrain Management>Add Feature Class to Terrain)
- If you don't have edge of water breaklines, you will want to delineate them yourself and add them to the terrain model. Add breaklines. (3D Analyst>Terrain Management>Add Feature Class to Terrain)
- Then build the terrain dataset. (3D Analyst>Terrain Management>Build Terrain)
After this step you should be able to view the terrain data on-screen and see if it is interpolating correctly over water features. If it doesn't look right, you may have to change the surface feature type for your breaklines. Here's the link to the help page.
Once you have everything looking right in the terrain dataset you can simply convert it to a DEM
(3D Analyst> Conversion> From Terrain> Terrain to Raster)
If you don't have 3D Analyst, you could digitize a region for your water body, convert it to a grid with a constant elevation, and then use the Con tool in Spatial Analyst to replace the values of the original DEM with the new values from your newly created water elevation grid.
- Digitize water polygon and specify the desired elevation in an attribute field.
- Convert to raster using the elevation field as the raster value. (Conversion Tools>To Raster>Polygon to Raster)
- Use the Con tool to replace the original values with the new water body elevation values. (Spatial Analyst Tools>Conditional>Con)