Imaging and prediction of fluid flow in the subsurface provides information that is crucial for decision making processes in fields such as groundwater management and enhanced oil recovery. The flow of an injected fluid through a reservoir depends primarily on the hydraulic conductivity, which is in general unknown or known only with low accuracy. A common way of imaging the flow is thus to intelligently modify the hydraulic conductivity model and simulate the fluid flow and geophysical imaging data that approximately match the observations over time. This process is also known as history matching. As the imaging process is a highly underdetermined inverse problem, we propose a new technique that avoids estimation of hydraulic conductivities. Instead, our approach directly estimates the flow field and initial distribution of the fluid from a time series of geophysical imaging data. Our method combines the flow equations with geophysical imaging to form a single inverse problem, where the unknowns are the initial state of the reservoir and the flow field. We discuss consistent discretization techniques, tailor specific regularizations, and use a modification of the variable projection method to solve the discrete optimization problem. We demonstrate the potential of our method on a model problem and show that our approach yields an improved flow estimate as well as an improved image quality. Finally, we show that the estimated flow field allows for the reconstruction of the subsurface structure.