Numerical solution of two-dimensional (2D) stochastic integral equations due to randomness has its own difficulties. For instance, most of them do not have analytical solution or obtaining their analytical solution is very hard. So, one of the essential requirement is presenting an efficient method to approximate the solutions of these equations with proper precision. In this paper, we introduce a numerical technique based on two dimensional modification of hat functions (2D MHFs) to estimate the numerical solution of 2D linear stochastic Volterra– Fredholm integral equations. In this approach, first the mentioned equation is transformed into a linear system of algebraic equations and then this system is solved via a direct or numerical method. Furthermore, it is proved that the order of convergence of the proposed method is O(h3). Finally, accuracy of this scheme is measured by solving two test examples via described technique.