Heterogeneous adsorbents, those composed of multiple surface and pore types, can result in stepwise isotherms that have been difficult to model. The complexity of these systems has often led to appealing to empirical equations without physical insights, unrealistic assumptions with many parameters, or applicability limited to a particular class of isotherms. Here, we present a statistical thermodynamic approach to model stepwise isotherms, those consisting of either an initial rise followed by a sigmoid or multiple sigmoidal steps, founded on the rigorous statistical thermodynamic theory of sorption. Our only postulates are (i) the finite ranged nature of the interface and (ii) the existence of several different types of microscopic interfacial subsystems that act independently in sorption. These two postulates have led to the superposition scheme of simple surface (i.e., Langmuir type) and cooperative isotherms. Our approach has successfully modeled the adsorption on micro-mesoporous carbons, gate-opening adsorbents, and hydrogen-bonded organic frameworks. In contrast to the previous models that start with a priori assumptions on sorption mechanisms, the advantages of our approach are that it can be applied universally under the above two postulates and that all of the fitting parameters can be interpreted with statistical thermodynamics, leading to clear insights on sorption mechanisms.