We consider the discretization of time-space diffusion equations with fractional derivatives in space and either 1D or 2D spatial domains. The use of implicit Euler scheme in time and finite differences or finite elements in space, leads to a sequence of dense large scale linear systems describing the behavior of the solution over a time interval. We prove that the coefficient matrices arising in the 1D context are rank structured and can be efficiently represented using hierarchical formats (H-matrices, HODLR). Quantitative estimates for the rank of the off-diagonal blocks of these matrices are presented. We analyze the use of HODLR arithmetic for solving the 1D case and we compare this strategy with existing methods that exploit the Toeplitz-like structure to precondition the GMRES iteration. The numerical tests demonstrate the convenience of the HODLR format when at least a reasonably low number of time steps is needed. Finally, we explain how these properties can be leveraged to design fast solvers for problems with 2D spatial domains that can be reformulated as matrix equations. The experiments show that the approach based on the