The ensemble Kalman filter (EnKF) is a powerful tool for assimilating data in earth system models. The approach allows real time Bayesian updating of system states and parameters as new data become available. This paper focuses on EnKF data assimilation in models of groundwater flow through complex geologic media. It has become common to treat the hydraulic conductivity of such media as correlated random fields conditioned on measured conductivity (medium property) and/or hydraulic head (system state) values. This renders the conductivity nonstationary and the corresponding conditional flow equations stochastic. Solving these equations and coupling them with EnKF generally entails computationally intensive Monte Carlo (MC) simulation. We propose to circumvent the need for MC through a direct solution of approximate nonlocal (integrodifferential) equations that govern the space-time evolution of conditional ensemble means (statistical expectations) and covariances of hydraulic heads and fluxes. We illustrate and explore our approach on synthetic two-dimensional examples in which a well pumps water from a randomly heterogeneous aquifer subject to prescribed head and flux boundary conditions. Embedding the solution in EnKF provides sequential updates of conductivity and head estimates throughout the space-time domain of interest. We demonstrate the computational feasibility and accuracy of our methodology, showing that hydraulic conductivity estimates are more sensitive to early than to later head values and improve with increasing assimilation frequency at early time.