Abstract One of the most serious limitations of current astrochemical models with the rate equation (RE) approach is that only a single type of binding site is considered in grain surface chemistry, although laboratory and quantum chemical studies have found that surfaces contain various binding sites with different potential energy depths. When various sites exist, adsorbed species can be trapped in deep potential sites, increasing the resident time on the surface. On the other hand, adsorbed species can be populated in shallow sites, activating thermal hopping and thus two-body reactions even at low temperatures, where the thermal hopping from deeper sites is not activated. Such behavior cannot be described by the conventional RE approach. In this work, I present a framework for incorporating various binding sites (i.e., binding energy distribution) in gas-ice astrochemical models as an extension of the conventional RE approach. I propose a simple method to estimate the probability density function (pdf) for the occupation of various sites by adsorbed species, assuming a quasi-steady state. By using thermal desorption and hopping rates weighted by the pdfs, the effect of binding energy distribution is incorporated into the RE approach without increasing the number of ordinary differential equations to be solved. This method is found to be accurate and computationally efficient, and enables us to consider binding energy distribution even for a large gas-ice chemical network which contains hundreds of icy species. The impact of the binding energy distribution on interstellar ice composition is discussed quantitatively for the first time.