As far as the knowledge of the seabed is concerned, both for safe navigation and for scientific research, 3D models, particularly digital bathymetric models (DBMs), are nowadays of fundamental importance. This work aimed to evaluate the quality of DBMs according to the interpolation methods applied to obtain grid format 3D surfaces from scattered sample points. Other complementary factors affecting DBM vertical accuracy, such as seabed morphological complexity and surveyed points sampling density, were also analyzed by using a factorial ANOVA experimental design. The experiments were performed on a multibeam dataset provided by the Italian Navy Hydrographic Institute (IIM) with an original resolution of 1 m × 1 m grid spacing, covering a surface of 0.24 km2. Six different sectors comprising different seabed morphologies were investigated. Eight sampling densities were randomly extracted from every sector, each with four repetitions. Finally, four different interpolation methods were tested, including: radial basis multiquadric function (RBMF), ordinary kriging (OK), universal kriging (UK) and Gaussian Markov random fields (GMRF). The results demonstrated that both RBMF and OK produced very accurate DBM in areas characterized by low levels of seabed ruggedness at sampling densities of only 0.0128 points/m2 (equivalent grid spacing of 8.84 m). In contrast, a higher density of 0.1024 points/m2 (3.13 m grid spacing) was required to produce accurate DBM in areas with more complex seabed topography. On the other hand, UK and GMRF were strongly influenced by morphology and sampling density, yielding higher vertical random errors and more prone to slightly overestimate seabed depths. In addition, sampling density and morphology were the factors that most influenced the vertical accuracy of the interpolated DBM. In this sense, the highly statistically significant influence of the interaction between sampling density and morphology on the vertical accuracy of the interpolated DBM confirms the need to perform a preliminary analysis of seabed morphological complexity in order to increase, if necessary, the number of surveyed points in cases of complex morphologies.