A geostatistical methodology is proposed for integrating elevation estimates derived from digital elevation models (DEMs) and elevation measurements of higher accuracy, e.g., elevation spot heights. The sparse elevation measurements (hard data) and the abundant DEM-reported elevations (soft data) are employed for modeling the unknown higher accuracy (reference) elevation surface in a way that properly reflects the relative reliability of the two sources of information. Stochastic conditional simulation is performed for generating alternative, equiprobable images (numerical models) of the unknown reference elevation surface using both hard and soft data. These numerical models reproduce the hard elevation data at their measurement locations, and a set of auto and crosscovariance models quantifying spatial correlation between data of the two sources of information at various spatial scales. From this set of alternative representations of the reference elevation, the probability that the unknown reference value is greater than that reported at each node in the DEM is determined. Joint uncertainty associated with spatial features observed in the DEM, e.g. the probability for an entire ridge existing, is also modeled from this set of alternative images. A case study illustrating the proposed conflation procedure is presented for a portion of a USGS one-degree DEM. It is suggested that maps of local probabilities for over or underestimation of the unknown reference elevation values from those reported in the DEM, and joint probability values attached to different spatial features, be provided to DEM users in addition to traditionally reported summary statistics used to quantify DEM accuracy. Such a metadata element would be a valuable tool for subsequent decision-making processes that are based on the DEM elevation surface, or for targeting areas where more accurate elevation measurements are required.