The accuracy of global tidal and mean sea surface (MSS) models degrades significantly in coastal and estuarine regions. These geophysical models are important for correcting measurements from satellite altimetry and are used in numerous scientific and engineering applications. The new Surface Water Ocean Topography (SWOT) mission is providing measurements at unprecedented horizontal resolution in these regions. This dataset provides both the opportunity and the need to quantify and correct the spatial variability of these errors. We develop a variational Bayesian framework for tidal harmonic analysis and MSS estimation which can be applied in the early stages of SWOT. The approach demonstrates superior robustness to different types of noise contamination in comparison to conventional least-squares approaches whilst providing full uncertainty estimation. By imposing a spatially coherent inductive bias on the model, we achieve superior harmonic constituent inference using extremely sparsely sampled data. Bayesian uncertainty estimation gives rise to statistical methods for outlier removal and constituent selection. Using our approach, we estimate a lower bound for the residual tidal variability for two SWOT Cal/Val passes (003 and 016) around the European Shelf to be at least 7% on average. We also show similar estimates cannot be produced using standard least-squares approaches. Tide gauge validation in the same region confirms the superiority of our empirical approach in coastal environments. Analysis of the estimated MSS error highlights the danger of interpolation and averaging in masking small-scale oceanographic processes. Empirical corrections for the SWOT data products are provided alongside an open-source Python package, VTide.