A conceptually simple general method to simulate model heterogeneous surfaces is proposed. It allows one to distribute adsorption energies over a lattice of adsorption sites according to an arbitrary energy distribution which can include spatial site energy correlations. The method provides a powerful tool to simulate either discrete or continuous site energy distributions with any kind of spatial correlations between site energies which makes it useful to test theoretical model of surface processes against MC simulations on the model surface. Effect of correlations on equilibrium and kinetics of gas adsorption thermal desorption, surface diffusion and reactions can be studied by MC simulation. More sophisticated model heterogeneous surfaces incorporating distributions of saddle point energies (consequently, activation barriers for diffusion) can also be treated. The method is the same as that used for, evaluating the equilibrium states of an ensemble of interacting particles through Monte Carlo simulation.