I want to find areas of relative "flatness" in a DEM.
How can I find flat areas?
I want to find a neighbourhood of contiguous cells that have an elevation variance that is within a certain tolerance threshold.
I can use the ESRI suite of tools, or any open source GIS.
If you're looking for a region of specific size and area:
This problem has many different solutions: it needs further restrictions. After all, by choosing any value within the range of elevations in the dataset, the contour for that value will be "flat" (and horizontal). This probably is not the sought-for answer, but it is an answer and it's perfectly valid. Moreover, among all possible answers, the contour lines will achieve minimal "elevation variance," so arguably they are the best answers--at least until additional criteria are specified.
Typically one restricts such searches to reasonably compact regions by asking that they be circular, rectangular, or perhaps by limiting their tortuosity. Also, there are many possible definitions of "elevation variance," such as range, standard deviation, IQR, MAD, etc.
As an example of what can easily be done, suppose we want to find circular regions where the elevation range does not exceed some threshold. (This is useful for implementing local requirements which frequently say a plot of land cannot rise more than z feet over distance x.) The calculation is dead simple: pick the radius you want to investigate and compute the focal range for a circular neighborhood of that radius. Every cell in the resulting grid whose value is within the elevation tolerance is the center of a "flat" circle.
There are obvious variations: you can explore a range of radii to search for the largest possible such circle; you can use neighborhoods of different shapes (such as square neighborhoods or even irregular neighborhoods); you can compute focal standard deviations instead of focal ranges. I believe all these options have been discussed at one time or another on the old ESRI forums.
One huge advantage this approach has over a slope-based solution is that it is immune to the small measurement errors that afflict almost all DEMs. Consider, for example, a DEM that was computed from a contour map. Unless the DEM processing was extremely well done, it will still show terrace-like artifacts near the former contour lines. At those cells, computed slopes can be extreme (but erroneous). This won't materially affect the elevation-based focal stats technique but it will cause slope-based techniques to give the wrong, or inferior, results.
Your best bet would be to first calculate slope, as om_henners has mentioned, then use a focal operator to collapse the result down to the scale you're interested in. I'd recommend using FocalMean, where the neighborhood size matches the scale of 'flat' areas you're interested in. So for example, if you had a raster of 10m resolution, a 5x5 neighborhood would give you the average slope over 2.5 km2 area.
The same analysis could be done in GRASS, with r.mfilter.
FocalStd directly on the DEM itself. Either way, what if I do not want to find only areas of a particular geometry and size? That is, I want areas that maximise the shape and size while fitting within the constraints.
As far as I understand your question, you might be very interested in the techniques mentioned in this paper: Gallant, J.C., Dowling, T.I. (2003): 'A multiresolution index of valley bottom flatness for mapping depositional areas', Water Resources Research, 39/12:1347-1359 and implemented in SAGA GIS: http://www.saga-gis.org/saga_modules_doc/ta_morphometry/ta_morphometry_08.html
If you have access to Spatial Analyst you could use Slope instead of cut/fill.
You could try Landserf.
There is a tool to compute slope maps:
(flat areas are in yellow)
See section 1.6 of the user guide.
An other possibility is to triangulate your DEM and compute the triangle's slopes.
