Surface-related volume calculations are usually performed on: (i) a single grid (which we call here the “upper” surface) with respect to a reference plane (i.e. z=constant, which we call the “lower” surface), rather as per the surface area computations described in Section 4.2.3, Surface area; or (ii) on a pair of grids. In the latter case the difference in volumes is computed assuming that the first grid represents the upper surface and the second represents the lower surface, even though the two may intersect. The result of the computations is an estimated positive volume (upper surface volume that is above the lower surface) and an estimated negative volume (where the lower surface is above the upper surface). If the upper surface is always above the lower then all values will be positive. Some packages require an additional reference surface to be defined, e.g. as an absolute base value.
For example, using our test surface TQ81NE, with a reference plane of z=0, the entire surface is above sea level so the total volume reported is entirely positive with a value of roughly 570 million cubic meters, i.e. around 0.57km3. Note that this figure includes adjustment to allow for a 10m x 10m grid size. Using a reference plane at z=30m the positive (solid matter) volume is around 0.13km3 and the negative (air space) around 0.3km3. At 22.5m the two volumes are approximately equalized at c.0.2km3, i.e. cut=fill if one were trying to level out this region using the available materials.
Volumes are calculated by numerical integration procedures, similar to those described earlier for determining the area of a polygon. With very small grid sizes the accuracy of volume estimation is, in theory, better, but this will depend on the way in which the surface is modeled. Surface roughness affects the results, especially where the surface is complex in form, the differences between upper and lower surface are small, and the computational method simplistic. The simplest procedure (which appears to be that adopted by ArcGIS) involves use of the trapezoidal rule (see also, Section 4.2.1, Length and area for vector data), first in the x-direction across all grid nodes or cells, for the first row in the grid, and then for each subsequent row. This creates a set of slices with areas A1,A2,…An. The slices times the grid resolution in the y-direction provides an estimate of volume. The first and last elements in each row may be weighted 1, with all other elements weighted 2, and then the sums divided by 2 to adjust for double counting. Although this provides a fairly crude and clearly biased estimate of the surface form at the top of each DEM cell, it is quite effective — the difference between this approach and more complex procedures, such as the extended Simpson’s rule and variants of this which are provided in Surfer for example, is minimal in the test example we have cited (a tiny fraction of a percent).
Engineering-oriented packages, such as Autodesk’s Land Desktop product, and civil engineering extensions to this, offer a wider range of options including volume estimation from surface point data and cross-sectional slices. ArcGIS includes volume computation from TIN surface representations, which one would presume is exact with respect to this surface model. Mathematical packages, such as MATLab perform surface integration in much the same way as described above, although the input is typically a mathematical function rather than a terrain model. This enables the numerical integration procedure to meet precise tolerances (accuracy levels) based on using successively finer steps (i.e. similar to utilizing an arbitrarily fine grid model).