Anomalous sub-diffusion processes governed by non-local operators have been used in various applications. These are recently adapted to model complex single-and multi-phase fluid flow in unconventional reservoirs to understand the long-time behavior of various quantities of interest. Such processes, requiring long-time simulation, are modeled by a non-local in-time fractional derivative based Darcy's law and a conservation of mass equation, leading to a coupled system of first-order fractional partial differential equations (FPDEs). Standard time-stepping discretization of these continuous models lead to computational bottleneck for long-time simulation. Memory limitations [in each node of high performance computing (HPC) clusters] dictate computational constraints for resolving fine scale spatial structures in the models. We describe and implement a parallel-in-time and parallel-in-space mixed finite element (FEM) based discretization computer model to simulate the FPDEs. Our hybrid parallel-in-time-and-space framework facilitates efficient long-time simulation of the first-order computer model to overcome the time-stepping computational bottleneck and, within the node memory constraints, approximate the continuous models with large degrees of freedom (DoF).