AbstractIn this article, we develop a Crank–Nicolson alternating direction implicit finite volume method for time‐dependent Riesz space‐fractional diffusion equation in two space dimensions. Norm‐based stability and convergence analysis are given to show that the developed method is unconditionally stable and of second‐order accuracy both in space and time. Furthermore, we develop a lossless matrix‐free fast conjugate gradient method for the implementation of the numerical scheme, which only has memory requirement and computational complexity per iteration with N being the total number of spatial unknowns. Several numerical experiments are presented to demonstrate the effectiveness and efficiency of the proposed scheme for large‐scale modeling and simulations.