Abstract. We present here an aerosol model selection based statistical method in Bayesian framework for retrieving atmospheric aerosol optical depth (AOD) and pixel-level uncertainty. Especially, we focus on to provide more realistic uncertainty estimate by taking into account a model selection problem when searching for the solution by fitting look-up table (LUT) models to a satellite measured top-of-atmosphere reflectance. By means of Bayesian model averaging over the best-fitting aerosol models we take into account an aerosol model selection uncertainty and get also a shared inference about AOD. Moreover, we acknowledge model discrepancy, i.e. forward model error, arising from approximations and assumptions done in forward model simulations. We have estimated the model discrepancy empirically by a statistical approach utilizing residuals of model fits. We use the measurements from the TROPOspheric Monitoring Instrument (TROPOMI) onboard the Sentinel-5 Precursor in ultraviolet and visible bands, and in one wavelength band 675 nm in near-infrared, in order to study the functioning of the retrieval in a broad wavelength range. We exploit a fundamental classification of the aerosol models as weakly absorbing, biomass burning and desert dust aerosols. For experimental purpose we have included some dust type of aerosols having non-spherical particle shapes. For this study we have created the aerosol model LUTs with radiative transfer simulations using the libRadtran software package. It is reasonably straightforward to experiment with different aerosol types and evaluate the most probable aerosol type by the model selection method. We demonstrate the method in wildfire and dust events in a pixel level. In addition, we have evaluated in detail the results against ground-based remote sensing data from the AErosol RObotic NETwork (AERONET). Based on the case studies the method has ability to identify the appropriate aerosol types, but in some wildfire cases the AOD is overestimated compared to the AERONET result. The resulting uncertainty when accounting for the model selection problem and the imperfect forward modelling is higher compared to uncertainty when only measurement error is included in an observation model, as can be expected.