I formulate a digital twin approach to the reconstruction of velocity fields from noisy and sparse magnetic resonance velocimetry signals. The method learns the most probable fluid dynamics model that fits the data by solving a Bayesian inverse Navier–Stokes boundary value problem. This jointly reconstructs and segments the velocity field, and at the same time infers hidden quantities such as the hydrodynamic pressure and the wall shear stress, as well as their uncertainties. Using a Bayesian framework, I regularize the problem by introducing a priori information about the unknown parameters in the form of Gaussian random fields. This prior information is updated using the Navier–Stokes problem, an energy-based segmentation functional, and by requiring that the reconstruction is consistent with the signals. I create an algorithm that solves this inverse problem, and first test it for noisy synthetic images of 2D flow in a simulated aortic aneurysm. I then extend the method to noisy, sparsely-sampled signals, and test it for experimental flow through a converging nozzle. I find that the method is capable of reconstructing and segmenting the velocity fields from sparsely-sampled (15% sampling), low (∼10) signal-to-noise ratio (SNR) signals, and that the reconstructed velocity field is almost identical to that derived from fully-sampled (100% sampling) high (>40) SNR signals of the same flow. Finally, I implement the same algorithm in 3D and test it for an experimental flow through a 3D-printed physical model of an aortic arch. I show that the method can successfully reconstruct noisy flows in realistic geometries and high Reynolds numbers. Further, the method naturally extends to 3D periodic and unsteady flows, and its computational complexity can be substantially decreased if an adaptive discretization method (e.g. wavelets) is used. The reconstruction of in vivo cardiovascular and porous media flows are among the most exciting applications of this work.