Two-dimensional Volterra-Fredholm integral equations of Hammerstein type are studied. Using the Banach Fixed Point Theorem, the existence and uniqueness of a solution to these equations in the space L ∞ ([0, 1] × [0, 1]) is proved. Then, the operational matrices of integration and product for two-variable Bernoulli polynomials are derived and utilized to reduce the solution of the considered problem to the solution of a system of nonlinear algebraic equations that can be solved by Newton's method. The error analysis is given and some examples are provided to illustrate the efficiency and accuracy of the method. Keywords Two-dimensional integral equations • Volterra-Fredholm integral equations of Hammerstein type • Bernoulli polynomials • Operational matrix method • Collocation method Mathematics Subject Classification 65R20 • 45D05 1 Introduction Two-dimensional (2D) integral equations appear in the mathematical modeling of various physical phenomena and engineering problems. For example, it is shown in McKee et al. (2000) that a type of nonlinear telegraph equations is equivalent to 2D Volterra integral equations (VIEs). Also, in Graham (1980/1981), an application of 2D Fredholm integral Communicated by Hui Liang.