This study focuses on the analysis of adsorption equilibria of phenol and caffeine using two chromatographic columns employed in hydrophilic interaction liquid chromatography (HILIC): Hypersil GOLD HILIC and Shodex SILICA 5NH 4D. The aim of the study was to determine the energetic properties of the surface of the tested adsorbents depending on the composition of the mobile phase and in systems with two chemically different analytes. Adsorption isotherms were determined experimentally in systems with mobile phases of varying proportions of acetonitrile and methanol. The obtained equilibrium data were first analyzed using Scatchard plots, and in the in the subsequent step, the adsorption energy distribution was determined using the EM (Expectation–Maximization) method. In both columns, low-energy adsorption sites play a dominant role. However, under specific chromatographic conditions, the involvement of high-energy adsorption sites is also observed, indicating the heterogeneous nature of the stationary phase surfaces. The results enabled the selection of the most appropriate adsorption isotherm model and confirm that its proper choice is crucial for the correct interpretation of adsorption equilibrium data in HILIC systems. At the same time, they highlight the importance of accounting for surface heterogeneity in the development and optimization of chromatographic methods.