The main aim of the present study is to expand the operational matrix method for solving the fractional two-dimensional nonlinear weakly singular partial integro-differential equations.  To do this, firstly, we use and present the operational matrix of fractional integration of two-dimensional Fibonacci polynomials.  Then, by using the obtained operational matrices to approximate the fractional derivative of the solution of the considered equation, we convert the original problem to a nonlinear system of equations. Also, we present the error analysis of the proposed method by a theorem. Finally, we present and solve some numerical examples to illustrate the proposed method.