Increased human demand on the marine environment and associated biodiversity threatens sustainable delivery of ecosystem goods and services, particularly for shallow shelf-sea habitats. As a result, more attention is being paid to quantifying the geographical range and distribution of seabed habitats and keystone species vulnerable to human pressures. In this study, we develop a workflow based on unsupervised K-Means classification units and Generalized Linear Models built from multi-frequency backscatter analyses (95, 300 kHz), bathymetry and bathymetry derivatives (slope) to predict different levels of sandeel densities in Hempton's Turbot Bank Special Area of Conservation (SAC). For Hyperoplus lanceolatus densities, the performance of single frequency verses multi-frequency models is compared. Relatively high agreement between K-Means clustering outputs (from 95 kHz and multi-frequency models) and ground-truthed sandeel densities is noted. Moreover, Root Mean Squared Error (RMSE) values in this instance demonstrate that single-frequency models are favoured over the multi-frequency model in terms of predictive ability. This is mostly linked to the species strong affinity for sedimentary environments whose variability is better captured by the lower frequency system. Generally, these results provide important information about species-habitat relationships and pinpoint bedform features where sandeels are likely to be found and whose variability is potentially linked to the bathymetry domain. The workflow developed in this study also provides a proof of concept to support the design of a robust species-specific monitoring plan in marine protected areas. Most importantly, we highlight how decisions made during sampling, data handling, analysis could impact the final outputs and interpretation of Species Distribution Models and benthic habitat mapping.