Fractional evolution equations can characterize the physical and engineering phenomena by simulating memory effects and non-locality [1,2,3,4,5]. The nonlinear form of equations of this type provides a more accurate description of certain reality phenomena in comparison with the linear form [6,7,8,9]. Many scholars have attempted to seek numerical solutions since analytical solutions of these equations are mostly challenging to obtain [10,11,12,13]. In this paper, 2D nonlinear time-dependent Volterra partial integro-differential equations with multi-term weakly singular kernels are considered as follows:

with the initial condition

and the Dirichlet boundary condition

where
$\mu >0$,
$\Delta$ represents the 2D Laplacian operator, the third-order spatiotemporal mixed partial derivative introduces the viscosity effect, and
$I^{\left(\alpha \right)}$ is the
$\alpha$-order Riemann–Liouville fractional integral operator defined as follows:

which is weakly singular concerning time, and
$\Gamma (·)$ is the Gamma function. If the integral items in Equation (1) are removed, then Equation (1) becomes a classical Sobolev-type equation, which has been the subject of study by many researchers [14,15].
$\Omega \subset {\mathbb{R}}^2$ is an open and bounded rectangular region with the corresponding boundary
$\partial \Omega$, the nonhomogeneous terms
$f(x,y,t)$,
$u_0(x,y)$, and
$u_b(x,y,t)$ are all known functions, and the nonlinear term
$G\left(u\right)\in C^2\left(\mathbb{R}\right)$ is Lipschitz continuous with the Lipschitz constant
$L>0$, that is,

Particularly, Equation (1) is linear when
$G\left(u\right):=u$ or 0. The linear situation of Equation (1) arises from many models for heat flow in a rectangular, orthotropic material with memory, which has been investigated by some researchers [16,17,18]. Recently, significant research has been carried out regarding problems analogous to (1). In [19], Chen et al. proposed a BDF2 compact difference scheme to solve a class of
$\alpha$-order Riemann–Liouville fractional integral equations and achieved the temporal convergence rate with
$1+\alpha$. In [20], Qiao et al. investigated a class of integro-differential equations with multi-term weakly singular kernels using the BDF2 ADI orthogonal spline collocation scheme, resulting in a convergence rate with
$\min 1+{\alpha }_1,1+{\alpha }_2$. Due to the insufficient smoothness of solution close to initial point, these methods do not succeed in attaining the expected second-order convergence, though their convergence orders are generally higher than first-order. The use of graded meshes with smaller time steps near
$t=0$ can handle the initial weak singularity of the solution [21]. Several other studies related to graded meshes are outlined in [22,23,24,25]. The product integration rule based on the piecewise linear interpolation is a second-order quadrature formula, which is commonly used to approximate integrals [26]. Since piecewise linear interpolation imposes no special requirements on the nodes, the graded meshes can be applied to the product integration rule. In [27], a BDF2 ADI OSC scheme on graded meshes was proposed to solve the 3D nonlinear fractional evolution equation, whose integral term is approximated by the product integration rule. Wu et al. in [28] constructed a second-order CN difference method on graded meshes to solve the fourth-order evolution equation with multi-term integrals. Their theoretical analyses referred to Lemma 4.2.3 in [19]; it is shown that the product integration rule achieves second-order convergence for suitable grading exponents. In this paper, we focus on resolving a general 2D nonlinear case of Equation (1) and handle the weak singularities of the solution at the initial moment based on the graded meshes. The smoothness assumptions of the solution follow the reference [29], that is,

where
$\alpha =\max {\alpha }_1,{\alpha }_2$, c refers to a generic positive constant, which is not necessarily uniform in distinct circumstances in this paper. Additionally, we perform a dimensional analysis on Equation (1). Suppose the dimensions of u, the time variable t, and the spatial variables x and y are represented by
$\left[U\right]$,
$\left[T\right]$, and
$\left[L\right]$, respectively, then

To ensure dimensional homogeneity, we assign the dimensions of the coefficients of the second, third, and fourth terms in Equation (1) as
${\left[L\right]}^2$,
${\displaystyle \frac{{\left[L\right]}^2}{{\left[T\right]}^{1+{\alpha }_1}}}$, and
${\displaystyle \frac{{\left[L\right]}^2}{{\left[T\right]}^{1+{\alpha }_2}}}$, respectively.

The non-negativity theorem presented by Lopez-Marcos in [29] is a forceful tool to testify to the stability and convergence of numerical approaches for fractional evolution equations; some specific applications can be found in [30,31,32]. In [33], Tang approximated the
$1/2$-order Riemann–Liouville fractional integral operator by the product integration rule and stated that the former
$n-1$ terms of the product integration rule satisfy the non-negativity theorem, while the n-th term does not. Furthermore, Chen et al. extended this result and proved that the first
$n-1$ terms of the product integration rule consistently satisfy the non-negativity theorem for any
$\alpha$ and graded meshes in [34]. However, numerical results demonstrate that it is a common phenomenon for the failure of the n-th term of the product integration method to satisfy the non-negativity theorem, which makes the proof complicated or unrigorous. In this paper, the stable and convergent properties of the numerical method are proved by modifying the n-th term and employing the discrete Gronwall lemma without affecting the whole algorithm.