In this paper, a new computational scheme based on operational matrices (OMs) of two-dimensional wavelets is proposed for the solution of variable-order (VO) fractional partial integro-differential equations (PIDEs). To accomplish this method, first OMs of integration and VO fractional derivative (FD) have been derived using two-dimensional Legendre wavelets. By implementing two-dimensional wavelets approximations and the OMs of integration and variable-order fractional derivative (VO-FD) along with collocation points, the VO fractional partial PIDEs are reduced into the system of algebraic equations. In addition to this, some useful theorems are discussed to establish the convergence analysis and error estimate of the proposed numerical technique. Furthermore, computational efficiency and applicability are examined through some illustrative examples. KEYWORDS two-dimensional Legendre wavelets, variable order fractional partial integro-differential equation, variable-order Caputo fractional derivative 1 INTRODUCTION In real life, many phenomena from diverse fields of science and engineering are modeled with nonlinear equations. The modeling of many dynamic systems leads to the integral term containing the unknown function when taking into account of memory effects. The integral and integro-differential