This paper analyzes a class of two-dimensional (2-D) time fractional reaction-subdiffusion equations with variable coefficients. The high-order L2-1 𝜎 time-stepping scheme on graded meshes is presented to deal with the weak singularity at the initial time t = 0, and the bilinear finite element method (FEM) on anisotropic meshes is used for spatial discretization. Using the modified discrete fractional Grönwall inequality, and combining the interpolation operator and the projection operator, the L 2 -norm error estimation and H 1 -norm superclose results are rigorously proved. The superconvergence result in the H 1 -norm is derived by applying the interpolation postprocessing technique. Finally, numerical examples are presented to verify the validation of our theoretical analysis.