Inspired by the experimental results by Haggerty, Zinn et al from 2004, we built an upscaled model for flow and coupled transport in heterogeneous media. The model extends those known from homogenization theory, and in particular those known as double porosity models, and applies to a wide range of scales and contrasts between media. The model includes three types of nonlocal (or memory) terms, whose kernels are derived from local cell problems. In the talk we discuss the challenges due to the presence of these memory terms and overview the recently proven strong, weak, and ultra weak stability of a numerical scheme. We also present comparison of the solutions to the microscale and macroscale models which demonstrate the relative importance of each of the three types of memory terms depending on the heterogeneity contrast. This is joint work with Ralph Showalter and Son-Young Yi.