Fractional reaction advection–diffusion problems have been of interest in the area of fractional calculus. In this paper, we consider two-dimension time and space nonlocal nonlinear reaction convection diffusion problem on a finite domain. We propose \(L_{1}\) interpolation approximation to approximate the Caputo fractional derivative combined the weighted and shifted-Grünwald Letnikov difference scheme to discretize the Riesz space fractional derivative. In addition, we use fully implicit Crank Nicolson method combined with weighted and shifted Grünwald–Letnikov difference operator and get a numerical approximation method is unconditionally stable and convergent with the accuracy of \({\mathcal {O}}(\tau ^{{2-\alpha }}+h_x^2+h_y^{2})\) . We take the alternative direction implicit approach for two-dimension problem, which is used to reduce a two-dimension problem to a one-dimension equation by applying x and y-direction alternately over the groundwater problem. Finally, the numerical simulations with examples of nonlinear source term are given, demonstrating the theoretical study and confirming the effectiveness of the proposed method.
