We test a fully non-linear method to solve Bayesian seismic tomographic problems using data consisting of observed travel times of first-arriving waves. Rather than using Monte Carlo methods to sample the posterior probability distribution that embodies the solution of the tomographic inverse problem, we use variational inference. Variational methods solve the Bayesian inference problem under an optimization framework by seeking the best approximation to the posterior distribution, while still providing fully probabilistic results. We introduce a new variational method for geophysics – normalizing flows. The method models the posterior distribution by employing a series of invertible and differentiable transforms – the flows. By optimizing the parameters of these transforms the flows are designed to convert a simple and analytically known probability distribution into a good approximation of the posterior distribution. Numerical examples show that normalizing flows can provide an accurate tomographic result including full uncertainty information while significantly decreasing the computational cost compared to Monte Carlo and other variational methods. In addition, this method provides analytic solutions for the posterior distribution rather than an ensemble of posterior samples. This opens the possibility that subsequent calculations that use the posterior distribution might be performed analytically.