Summary: In this paper, we focus on the numerical simulation for two-dimensional nonlinear fractional sub-diffusion equations in the presence of time delay. Firstly, we investigate the existence, uniqueness and regularity of the solution for such problems. The theoretical result implies that the solution at \(\tau^+\) is smoother than that at \(0^+\), where \(\tau\) is a constant time delay, and this is an improvement for the work [T. Tan et al., J. Sci. Comput. 92, No. 3, Paper No. 98, 26 p. (2022; Zbl 07568990)]. Secondly, a high-order difference scheme based on \(L1\) method is constructed. For the sake of repairing the convergence order in temporal direction and improving the computational efficiency, an efficient time two-grid algorithm based on nonuniform meshes is first developed. The convergence order of the two-grid scheme reaches \(O(N_F^{-\min\{r\alpha,2-\alpha\}}+N_C^{-\min\{2r\alpha,4-2\alpha\}}+h_1^2+h_2^2)\), where \(N_F\) and \(N_C\) represent the number of the fine and coarse grids respectively, while \(h_1\) and \(h_2\) are the space-step sizes. Furthermore, stability and convergence analysis of the proposed scheme are carefully verified by energy method. Finally, numerical experiments are carried out to show the validity of theoretical statements.