Determination of interfacial areas is crucial for accurate identification of multiphase subsurface flow processes, e.g., in situ water treatment, contaminant transport, phase change, convection/diffusion, and colloid adsorption/desorption/Migration. Pore-scale imaging and simulation techniques provide an appealing opportunity to explicitly predict such geometrical, characteristics of saturated and unsaturated porous media. Here, synthetic unconsolidated sand packs in a wide range of grain. arrangements and size distributions ate reconstructed to,extensively capture, all possible heterogeneous pore-level geometries, in the McMurray Formation, the primary bitumen formation in the Athabasca oil sands deposit. Multiphase fluid occupancies throughout the partially saturated media during drainage and imbibition are predicted applying a direct pixel-wised pore morphological model, incorporating a wetting phase layer covering rock/solid surfaces. The postprocessing results are verified using two-phase experimental data points and images, demonstratitig a remarkable variation of the interfacial area,as a function of saturation profile, rock-configuration, and displacement scenarios. Empirical models are proposed to predict the bulk and meniscus areas using average particle diameter and porosity as input parameters.