We develop a statistical framework for conducting inference on collections of time-varying covariance operators (covariance flows) over a general, possibly infinite dimensional, Hilbert space. We model the intrinsically non-linear structure of covariances by means of the Bures-Wasserstein metric geometry. We make use of the Riemmanian-like structure induced by this metric to define a notion of mean and covariance of a random flow, and develop an associated Karhunen-Loève expansion. We then treat the problem of estimation and construction of functional principal components from a finite collection of covariance flows, observed fully or irregularly. Our theoretical results are motivated by modern problems in functional data analysis, where one observes operator-valued random processes -- for instance when analysing dynamic functional connectivity and fMRI data, or when analysing multiple functional time series in the frequency domain. Nevertheless, our framework is also novel in the finite-dimensions (matrix case), and we demonstrate what simplifications can be afforded then. We illustrate our methodology by means of simulations and data analyses.