Quantitative mapping of seafloor sediment properties (eg. grain size) requires the input of comprehensive Multi-Beam Echo Sounder (MBES) datasets along with adequate ground truth for establishing a functional relation between them. MBES surveys in extensive shallow shelf areas can be a rather challenging and time-consuming task, resulting in time and cost intensive data collection. It is therefore often the problem of dealing with sparse data and/or data without full coverage. This study deals with MBES data acquired within an area of approximately 30 x 55 km (Cleaverbank, North Sea, the Netherlands) with a line spacing of ~2 km in 50 m average water depth. Additionally ground truth samples were taken with a Hamon grab in a regular grid of ~1 x 1 km with samples at and in-between some of the MBES survey lines. These ground truth data cover a subset of the area under investigation from variable depth and sediment types ranging from silty clay to boulders. Here we combine geostatistical and multivariate modelling for predictive mapping of the median grain size across the whole area. First bathymetric data was cleaned and raw backscatter values were classified into 7 classes representing different sediment types using a Bayesian method. A multiple linear regression was performed with the ground truth (sediment grain size) data resulting in an empirical function that shows high correlation (90%) of the median grain size with depth and backscatter classes. However missing information in un-surveyed areas needs to be estimated before predictive mapping can take place. For this purpose we interpolated sparse bathymetric and backscatter class data in order to obtain layers of continuous information for the whole area. To achieve this, MBES lines were treated as Single Beam Echo-Sounder data by producing a reduced, equally-spaced grid of bathymetry and backscatter class points, for each MBES line. Interpolation of this bathymetric data set using Ordinary Kriging (OK) yielded satisfactory results with low root-mean-square errors. For the backscatter class we utilized both OK and a Bayesian Maximum Entropy (BME) algorithm which is regarded as more suitable for categorical data. As backscatter classes are ordered by increasing acoustic intensity makes the interpolation results useful as they also include values between classes. Interpolated bathymetry and backscatter classes were validated with additional data from crossing lines. Subsequently, the interpolated layers produced by geostatistical modelling became input to the empirical regression function connecting median grain size and MBES data. The predictive layer of median grain size was assessed by using independent ground truth data from not-mapped locations. By applying linear regression to the interpolated data it was found that both, OK and BME provided results that predict median grain sizes with a high correlation (up to r2=0.85) to the observed values. The overall approach suggests that the combination of geostatistical with multivariate modelling can offer good results for large scale seafloor classification when data are limited. Additionally the prediction accuracy for each layer (bathymetry, backscatter class) can be estimated and spurious results can be identified. @en