This paper considers Bayesian inversion of seabed reflection-coefficient data for multi-layer geoacoustic models at several sites, with the goal of studying lateral variability of the seabed. Rigorous uncertainty estimation is carried out to resolve lateral variability of the sediments from inherent inversion uncertainties. The uncertainty analysis includes Bayesian model selection, comprehensive quantification of data error statistics, and a Markov-chain Monte Carlo approach to transforming data uncertainties to model uncertainties. Model selection is addressed using the Bayesian information criterion to ensure parsimony of the parametrizations. Data error statistics are quantified by estimating full covariance matrices from data residuals, with posterior statistical validation. A Metropolis-Hastings sampling algorithm is used to compute posterior probability densities. Four experiment sites are considered along a track located on the Malta Plateau, Mediterranean Sea, and the inversion results are compared to cores taken at each site. Differences between profile marginal-probability distributions at adjacent sites are quantified using the Bhattacharyya coefficient. Differences that exceed the estimated geoacoustic uncertainties are interpreted as spatial variability of the seabed. The results are compared to an interpretation of geologic features evident in a chirp sub-bottom-profiler section.