In fully water-saturated, porous geomaterials, such as fine-grained soils and cement mortars with fine pore spaces, freezing leads to the formation of ice within the pore network and the development of localized ice lenses. With a continued water supply, these ice lenses grow, inducing internal stresses that can damage the geomaterial skeleton. Therefore, understanding the initiation and propagation of ice lenses is of significant engineering interest for mitigating frost-induced damage in geomaterials. In this contribution, we present a computational, phase-field, fracture-based, thermo-hydro-mechanical (THM) finite element model to predict ice lens formation in geomaterials during freezing. Ice lens formation is modeled using a phase-field fracture formulation following the framework proposed in \cite{Heider2022,Suh2022}. Within this framework, we examine several theoretical criteria for triggering ice lens initiation as proposed in the literature \cite{Sweidan2022}. We use this model to investigate ice lens formation by numerically analyzing mesoscale soil samples. We explicitly account for porosity heterogeneity via random field representations, as presented in \cite{Andrade2008}. This study's central hypothesis is that the combined effect of porosity-driven material heterogeneity and an appropriate ice lens initiation criterion enables more realistic modeling of complex ice lens formation in geomaterials. To evaluate the predictive ability of our approach, we compare our numerical results with experimental observations reported in the literature for fine-grained soils.