SUMMARY Estimating subsurface properties from geophysical measurements is a common inverse problem. Several Bayesian methods currently aim to find the solution to a geophysical inverse problem and quantify its uncertainty. However, most geophysical applications exhibit more than one plausible solution. Here, we propose a multimodal variational autoencoder model that employs a mixture of truncated Gaussian densities to provide multiple solutions, along with their probability of occurrence and a quantification of their uncertainty. This autoencoder is assembled with an encoder and a decoder, where the first one provides a mixture of truncated Gaussian densities from a neural network, and the second is the numerical solution of the forward problem given by the geophysical approach. The proposed method is illustrated with a 1-D magnetotelluric inverse problem and recovers multiple plausible solutions with different uncertainty quantification maps and probabilities that are in agreement with known physical observations.