Geophysical inversion can produce 3D models of the subsurface's physical properties. The smoothness of property variations in these models makes it challenging to automatically find boundaries of homogeneous regions, where these boundaries may have implications for petrophysical transition and are significant for geologic interpretation. We have developed a new boundary detection technique that nonparametrically identifies and subtracts homogeneous regions from the 3D model, leaving exposed edges. The method is based on kernel density estimation of local property variations, in which the number of modes in the kernel density estimates (i.e., the local mode cardinality [MC]) is used to identify edge voxels within the model. Two edge detection operators were developed: one using local values exclusively and the other incorporating the spatial distribution of local values into the local MC, which ismore sensitive to local variations. To assist in the geologic interpretation, continuous boundary surfaces are generated from the identified edge voxels using a gradient-based 3D image morphological operation. The technique was evaluated on synthetic rock property models and effectively identified the lithologic boundaries, even with non-Gaussian noise. The proposed operators were also applied to visualize edges in density, conductivity, magnetic susceptibility, and seismic tomography models of the Kevitsa Ni-Cu-PGE deposit (Lapland, Finland) generated by geophysical inversion. The boundaries detected via the proposed technique can be used for visualization and may be useful in further geostatistical computation due to their statistical foundation.