In this work, we design a multiscale simulation method based on the Generalized Multiscale Finite Element Method (GMsFEM) for numerical modeling of fluid seepage under permafrost condition in heterogeneous soils. The complex multiphysical model consists of the coupled Richards equation and the Stefan problem. These problems often contain heterogeneities due to variations of soil properties. For this reason, we design coarse-grid spaces for the multiphysical problem and design special algorithms for solving the overall problem. A numerical method has been tested on two- and three-dimensional model problems. A a quasi-real geometry with a complex surface is considered for the three-dimensional case. We demonstrate the efficiency and accuracy of the proposed method using several representative numerical results.