Has probado rasterizar los vectores con valores de cero y uno (cero para los lagos, el interior) y luego multiplicar por el raster del MDE? De esta manera todos los valores al interior del lago quedarian con el valor de cero, mientras los otros valores se mantendrian originales, y se entenderian como una superficie plana, sin embargo los artefactos de la red de drenaje seguiran apareciendo en los alrededores del lago. Adicionalmente, podrias modificar el valor de NODATA, o incluir este en el calculo, me parece que en ArcGIS usan -9999 para NODATA, para que no sean considerados en el calculo.