  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 for the occupation of various sites by adsorbed
species, assuming a quasi-steady state. By using thermal desorption and hopping
rates weighted by the probability density functions, 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.
