## 1 Introduction  
At high altitudes and latitudes, the rock interior is rich in crystalline ice [1]. As human activity in the region increases, so do the demands on the structural properties of frozen rock [2; 3; 4]. In order to safely and efficiently carry out engineering works in permafrost regions, it is essential to study the mechanical properties and failure modes of frozen rocks. However, there have been few studies on the stress field distribution of frozen rocks, and methods for calculating the stress intensity factor of cracks in frozen rocks under biaxial compression are lacking. Therefore, it is crucial to construct a realistic ice-rock model and study the distribution of stress fields within the model, as well as methods for calibrating fracture parameters. In the field of mechanical modeling, layered structures have become the preferred choice for scientists due to their ability to provide uniform stress distribution. In layered structures, scholars considered the presence of ice as a nodular formation with a certain thickness and conducted a series of experimental studies. For example, [PERSON] [5] studied the variability of the thickness of intermediate rock layers in response to earthquakes,while [PERSON] [6] investigated the characteristics of large deformation in layered soft rock and proposed corresponding control techniques. In the field of simulation research, [PERSON] [7] utilized FLAC3D to conduct a multi-layer-rock true triaxial experiment. They investigated the anisotropy under real stratigraphic conditions and proposed the type of damage influenced by the rock matrix and weak planes. [PERSON] [8] analyzed the impact of formation thickness on the fracturing process of reservoirs using FD numerical simulation. These studies advanced the research progress regarding the mechanical properties of ice-rock. However, finding uniformly layered and continuous ice-rock samples in practical engineering is challenging, which imposes significant limitations on the current research model. In terms of the theoretical study of fracture parameters, scholars' research began with the analysis of an infinitely large flat plate subjected to uniform tensile stress. They have developed fruitful calculation methods, such as [PERSON]'s function method [9; 10] and [PERSON]'s series expansion method [11; 12]. Subsequently, the interaction between multiple holes and multiple cracks gradually became a popular research field. For example, [PERSON] [13] calculated the artificial loads on the surface of the cracks and circular holes by introducing the stress intensity factor to account for the interaction. He then derived a system of Cauchy-type integral equations. Finally, various methods, such as the displacement discontinuity method [14], the dislocation density method [15], and the boundary integral method [16], were introduced. For example, [PERSON] et al. [17] proposed a formula that utilizes the superposition principle to calculate the stress intensity factor of multiple holes and cracks under the influence of far-field stresses and surface stresses. [PERSON] [18] calculated the stress intensity factor at the tip of a closed-type co-linear crack but did not explain the significance of the covariance and did not provide further analysis of the stress and displacement fields at the crack tip. Due to the compressive loading in the far field, there is contact between the surfaces that are cracked. This causes some difficulties in stress analysis at the crack tip. In this regard, scholars gradually adopted computer-aided engineering (CAE) methods for further research. For example, [PERSON] [19] and [PERSON] [20] utilized different iterative algorithms to investigate the issue of contact friction. Additionally, [PERSON] [21] employed the extended finite element method in conjunction with augmented Lagrange multipliers to simulate the expansion process of contact friction cracks. This provides a basis for stress analysis and predicting the extension of closed-type cracks.