This research proposes a hybrid numerical method, for modeling the development of artificial ground freezing in soil containing subsurface inhomogeneities. In this approach, by combining the Boundary Element Method (BEM) and employing time-independent fundamental solutions, appropriate boundary integral equations are developed. The volume integrals arising from the presence of dynamic terms are incorporated into the equations using the Finite Element Method (FEM). For this purpose, a type of boundary-finite element, combining a quadratic boundary element with a three-node triangular finite element, was developed, and the solvable forms of the final equations were presented.

Subsequently, by implementing the hybrid method into a computational algorithm, its accuracy and efficiency was evaluated and validated by solving several benchmark examples. Finally, in a parametric study, the application of the hybrid method for modeling the development of artificial freezing in saturated soil containing a inhomogeneities in the form of an unsaturated rock mass is described. The effects of varying the cross-sectional area of the inhomogeneities and its distance from the freeze pipe were evaluated.

The results of the parametric study indicate that the presence of the inhomogeneities reduces the volume of the freeze bulb by up to 20%. Furthermore, inhomogeneities with a circular cross-section were more limited the development of freezing compared to a square shape.