Full waveform inversion (FWI) creates high resolution models of the Earth’s subsurface structures from seismic waveform data. Due to the non-linearity and non-uniqueness of FWI problems, finding globally best-fitting model solutions is not necessarily desirable since they fit noise as well as signal in the data. Bayesian FWI calculates a so-called posterior probability distribution function, which describes all possible model solutions and their uncertainties. In this paper, we solve Bayesian FWI using variational inference and propose a new methodology called physically structured variational inference, in which a physics-based structure is imposed on the variational distribution. In a simple example motivated by prior information from past FWI solutions, we include parameter correlations between pairs of spatial locations within a dominant wavelength of each other, and set other correlations to zero. This makes the method far more efficient in terms of both memory requirements and computation, at the cost of some loss of generality in the solution found. We demonstrate the proposed method with a 2D acoustic FWI scenario, and compare the results with those obtained using other methods. This verifies that the method can produce accurate statistical information about the posterior distribution with hugely improved efficiency (in our FWI example, 1 order of magnitude in computation). We further demonstrate that despite the possible reduction in generality of the solution, the posterior uncertainties can be used to solve post-inversion interrogation problems connected to estimating volumes of subsurface reservoirs and of stored CO2, with minimal bias, creating a highly efficient FWI-based decision-making workflow.