We consider a Bayesian multi-model fault slip estimation (BMMFSE), in which many candidates of the underground-structure model characterized by a prior probability density function (PDF) are retained for a fully Bayesian estimation of fault slip distribution to manage model uncertainty properly. We performed geodetic data inversions to estimate slip distribution in long-term slow slip events (L-SSEs) that occurred beneath the Bungo Channel, southwest Japan, in around 2010 and 2018, focusing on the two advantages of BMMFSE: First, it allows for estimating slip distribution without introducing strong prior information such as smoothing constraints, handling an ill-posed inverse problem by combining a full Bayesian inference and accurate consideration of model uncertainty to avoid overfitting; second, the posterior PDF for the underground structure is also obtained through a fault slip estimation, which enables the estimation of sequential events by reducing the model uncertainty. The estimated slip distribution obtained using BMMFSE agreed better with the distribution of deep tectonic tremors at the down-dip side of the main rupture area than those obtained based on strong prior constraints in terms of the spatial distribution of the Coulomb failure stress change. This finding suggests a mechanical relationship between the L-SSE and the synchronized tremors. The use of the posterior PDF for the underground structure updated by the estimation for the 2010 L-SSE as an input of the analysis for the one in 2018 resulted in a more preferable Bayesian inference, indicated by a smaller value of an information criterion.