Summarizing flatness of county with statistic Elevation variance per square mile

You might want to check out the concept of rugosity, specifically the Arc-chord ratio (Du Preez, 2015). Basically, you take the area of a triangulation of the surface and the plane-of-best-fit to the edges of the surface.

Note that the fractal dimension is important: as the surface resolution increases, the surface area of its triangulation approaches infinity, so your entire dataset needs to have the same resolution. You don't mention which dataset you're using but SRTM would be ideal. The solution is robust to instances, like that mentioned by a commenter, where there's a single, tall structure on a mostly-flat terrain.

The algorithm should be fairly simple to implement in R (though I'm not able to post example code -- I've implemented it in C++ using the same general approach):

1. Convert the cells in the county to points, with the elevation as Z.
2. Perform a 3D Delaunay triangulation on the points.
3. Sum the area of the facets of the triangulation (one side only, obvs.)
4. Extract the points from the exterior ring of the triangulation and compute the plane-of-best-fit for those points, with the same bounding polygon as the triangulation.
5. Take the ratio of the areas of the triangulated surface and the POBF. This is the ACR.

One issue I foresee is the the memory requirement for geometric operations on the county level. You might have to down-sample the elevation raster or break up the work into smaller chunks and aggregate the result. You could run the entire country at some sampling resolution and then perform zonal statistics on the result by county.

(Note: I'm not certain that this is a good strategy and look forward to hearing from other members why that might or might not be so.)