Numerical investigation on origin and evolution of polygonal cracks on rock surfaces

We studied the formation and evolution mechanism of polygonal cracks on rock surfaces under cooling by modelling meso-damage mechanics, continuum mechanics and thermodynamics. Factors that affect rock surface damage include ambient temperature, lithology difference and boundary restrictions. We established and simulated a heterogeneous model with a surface weak layer for three types of boundaries. These were biaxial constraint, uniaxial constraint and free boundary. The initiation and propagation of polygonal cracks were reproduced for varying thickness and homogeneity of the weak layer. The results show that the boundary constraints strongly influence the polygonal cracking. Many polygonal or parallel cracks are generated on the rock surface under biaxial or uniaxial constraint. The unconstrained rock surface displays polygonal cracks at the center and parallel cracks in the surrounding areas. The thicker the surface weak layer, the larger the average area of formed blocks. Small blocks and short cracks are more numerous than large blocks and long cracks. As the heterogeneity index increases, the rock layer is more likely to produce blocks with relatively regular shapes. Quadrilateral, pentagonal and hexagonal blocks dominate regardless of changes in layer thickness and heterogeneity. However, the number of edges of the polygonal blocks is sensitive to rock heterogeneity. The polygons tend to become more complex with increasing inhomogeneity. This study contributes to understanding the complex formation mechanisms of polygonal cracks on rock surfaces in nature. Additionally, the simulations of three-dimensional fracture geometry provide a basis for developing algorithms to generate discrete fractures and blocks in discrete fracture network (DFN) analyses.


Introduction
Polygonal crack patterns are often observed on rock surfaces in nature. Most of these patterns are quadrilateral, pentagonal or hexagonal. Rock surfaces with well-developed cracks resemble tortoise shells, and are therefore termed 'craquelures' (Fig. 1). Rock craquelures occur widely on Earth. However, there are still disputes about their formation mechanism. When a rock surface experiences cooling-or drying shrinkage, high stresses concentrate at micro-defects and induce polygonal crack formation under certain circumstances. Similar mechanisms are responsible for the cracking of ceramics, oil paintings, murals and basalt columns, etc.
Rock is one of the most common materials used to construct the large-scale infrastructures that underpin human production and life. The safety and stability of engineering projects such as mines, tunnels, dams, and subgrades are significantly affected by the deformation and failure of rock masses. Thus, it is crucial to understand the formation mechanism of polygonal cracks and to understand the complex physical and mechanical behavior of rock masses for disaster prevention and mitigation. The cracking of rock surfaces is greatly influenced by environment, rock properties and mechanical state. Stress variations caused by changes in environment temperature and humidity can lead to rock fracturing (Schulke, 1973;Bradley et al., 1978). The effect of regional expansion and contraction has been investigated by many researchers (Jonhson, 1927;Rice, 1976;Smith, 1977;Williams and Robinson, 1989;Hall and André, 2003;Gómez-Heras et al., 2008;Riley et al., 2012;García-Rodríguez et al., 2015a). Siegesmund et al. (2000Siegesmund et al. ( , 2021 carried out the pioneering work on marble cracking caused by thermal expansion. Meanwhile, the thermal processes at environmental temperatures differ from those at higher temperatures because different temperature variations lead to diverse stress distributions due to mineral thermal expansion/contraction (Vázquez et al., 2015). Besides, the effect of the geometric shape of polygonal cracks and stress state cannot be ignored (Robinson and Williams, 1987;Riley et al., 2012;García-Rodríguez et al., 2015a). Considering that the cooling phases may be more effective in leading to stress redistribution which could result in crack propagation than heating phases due to the exponential change according to the Newton's cooling law (Gómez-Heras et al., 2008), this study mainly focuses on surface cooling.
On the surface of a vertical rock wall, the development of polygonal cracks is sometimes upward. Cracks intersect in a certain angle and may form square or rectangular polygons (Robinson and Williams, 1987). This orthogonal relation can be explained by a uniaxial stress state (Lachenbrunch, 1966). Tang et al. (2006) found that when the loading changed from uniaxial to biaxial, the parallel cracks change to polygonal cracks. They argued that equidistant cracks can be regarded as a special form of polygonal cracks. In addition, the formation of polygonal cracks is closely related to heterogeneity of the rock (Lachenbrunch, 1966;Hornig et al., 1996).
At present, most studies are based on observations of existing rock cracks. There is still a lack of systematic research on the geometric characteristics, fracturing modes and mechanical causes.
As regards theoretical research, some researchers (Zapperi et al., 1997;Tang et al., 1998;Liang, 2005) applied damage mechanics to study the damage mechanism from the micro-characteristics of rocks, to establish corresponding analysis models, and to further extend the relevant conclusions to the general brittle damage problem. Now, meso-damage mechanics has become a hot topic of damage theory. Some researchers (Zhu, 2000;Zhu and Tang, 2002;Chen et al., 2022;Tang, 2009;Feng et al., 2022; randomly assigned the mechanical properties of rocks according to a given statistical distribution and established random mechanical property models using the finite element method (FEM). However, the meso-mechanism of polygonal cracking is still not fully understood. Some research has been done to reveal crack initiation and development process. Morris et al. (1992) believed that the existence of microfractures in homogeneous materials made the tensile stress of materials concentrate at the tips, induced cracks to expand. The fracturing process can be considered as the expansion and connection of these microfractures under the action of stress. Polygonal cracking is closely related to the shrinkage of rocks. It is also affected by the mineral type and mineral content (Mitchell and Soga, 1976;Towner, 1987;Morris et al., 1994;Park et al., 2017). The properties of geomaterials, such as Poisson's ratio, the elastic modulus and the strength, have a great influence on polygonal cracking. Weinberger (1999) suggested that fractures might form below the surface first, and then expand upward until they intersect with surrounding fractures. Tang et al. (2008Tang et al. ( , 2010 and Wang et al. (2018) carried out a series of laboratory tests on ceramic crazing, and proposed a relationship between the occurrence and development of crazing and influencing factors including temperature, thickness, clay content and dry-wet cycle times. Simultaneously, they developed the Crack Image Analysis System (CIAS) method for analyzing the geometric characteristics of fracture networks (Tang et al., 2008). The CIAS (http://www.climate-engeo.com) is also used in this paper for image processing. However, due to the limitation of physical experiments, the gradual evolution of the stress fields leading to crazing was not discussed. Bao (2018) pointed out that the stress state between adjacent fractures experienced continuous evolution from tensile-to compressive stress, and that the transformation of stress state was the key reason for fracture saturation. However, the study of Bao (2018) was based on a two-dimensional simulation which ignored the spatial morphology of cracks and the interaction between cracks in threedimensional space.
The structure and distribution of polygonal cracks governs the deformation and strength of rock masses. Also, the geological history and stress state of rock masses can be determined by quantifying the crack features on a rock surface. Meanwhile, the complex formation mechanism of polygonal cracks can be revealed by understanding the core scientific problems of crack insertion and saturation. It is challenging to clarify the origin and evolution mechanism of polygonal cracks on rock surfaces (Tang et al., 2016;Wang and Konietzky, 2019;Saksala, 2019;Pressacco and Saksala, 2020).
In this paper, by combining meso-damage mechanics, continuum mechanics and thermodynamics, a heterogeneous model reflecting the characteristics of surface weak layers has been established. The formation and propagation of polygonal cracks on rock surfaces was reproduced numerically under three types of boundaries including biaxial constraint, uniaxial constraint and free boundary. The structure and distribution of polygonal cracking was investigated quantitatively for varying thicknesses and homogeneities of the weak layer. Besides, we discuss the morphology characteristics and mechanical mechanism of surface cracking. In addition, the value of quantifying the controls on three-dimensional (3D) fracture patterns and block geometry is for e.g. kinematic analysis of rock slope and underground stability, and to define fracture geometry for discrete fracture network (DFN) applications.

Constitutive relation
We used the 3D rock fracture process analysis (RFPA) method. To take rock heterogeneity into account, the relevant rock materials parameters were assigned according to the Weibull distribution function shown in Eq. (1).
where f(u) is the statistic quantity of the material parameter u; u is the element parameter (such as strength, elastic modulus, thermal expansion coefficient, etc.); u 0 is the average value of the parameter u; and m is the shape function of the distribution function, which can be termed the heterogeneity index. Besides, the heterogeneity index is related to rock mineralogy, and varies from 1 to the positive infinity depending on mineral type and mineral content. The higher the heterogeneity index m, the more uniform the material is. It should be determined by experiment for a specific kind of rocks.
In the RFPA method, the stress-strain relationship of a mesoscopic element under uniaxial state is as shown in Fig. 2. When the stress state of the element satisfies the specific failure criteria, the element will become damaged. The cumulative damage expression of the elastic modulus is: where E is the value of the elastic modulus after damage; E 0 is the initial value of the elastic modulus; and ω is the damage variable.
When the tensile stress borne by one mesoscopic element meets the corresponding tensile strength criterion, the damage variable ω is defined as: where λ is the residual strength coefficient which is defined as λ = f tr /f t0 ; ε t0 is the tensile strain corresponding to the ultimate elastic state and the initial damage threshold of the element; and ε tu is the ultimate tensile strain which is also related to the complete damage threshold of the element. When the compressional-shear stress sustained by the mesoscopic element satisfies the Mohr-Coulomb criterion, compressional-shear damage occurs and the damage variable ω can be determined as: where ε c0 is the ultimate compressive strain.

Temperature equation
The general thermodynamic coupling equations are: In these equations, σ ij and ε ij are stress and strain terms; F bi is mass force; ü i is the inertial variable; β is thermal stress coefficient; ΔT is the temperatures difference, ΔT = T − T 0 ; δ ij is Kronecker function; Q，λ， G，α，k，ρ and c are heat generation, Lame constant, shear modulus, thermal expansion coefficient, thermal conductivity, density and specific heat, respectively. These equations are interrelated and are solved together. However, for most engineering problems, the effect of strain on the temperature distribution is relatively small, and the inertial variable terms are also very small. That is, ε kk and ü i can be ignored. In this way, calculation of the temperature field only belongs to the heat conduction range, and calculation of the stress field includes the thermal stress. As a result, the temperature field can be calculated first. After that the stress and deformation caused by temperature and external load can be calculated.
Based on the above analysis, the thermal-mechanical coupling in this numerical model can be solved in two steps. First, the temperature field in the model is calculated according to the temperature conditions and second, the stress field is solved according to the mechanical conditions and the corresponding failure treatment is carried out. That is, each unit satisfies the following heat conduction equation and stress field equation.
The heat conduction equations include: The stress-strain field equations include: where k x , k y , k z , are the heat conduction coefficients of the unit in the x, y, and z directions, respectively; S 1 e , S 2 e , S 3 e are the 1st, 2nd, and 3rd thermal boundary conditions, respectively; t 0 e is the initial moment; φ(P, t)is the distribution function of temperature on the boundary point P with time t; q S2 is the boundary heat flow density; h is the heat transfer coefficient between the rock boundary and the outside world; T s is the ambient temperature; Ω e is the entire solution domain of the element; f i is the stress boundary; ũ i is the displacement boundary. The above equations are defined in the linear elastic range, but the macroscopic effects of rock deformation are nonlinear. How to handle the nonlinear characteristics of rock macroscopic deformation has long been studied.
In much previous work, most of the nonlinear deformation and fracture process of rock under stress has been attributed to elastic-plastic theory, which is expressed by macro-elastic-plastic theory. From a mesoscopic perspective, the nonlinearity of rock deformation is caused by the heterogeneity of the rock medium. The nonlinear nature of rock deformation is caused by continuous failure inside in rock during loading. Due to the extreme inhomogeneity of rocks, their properties vary greatly from the macro-to the micro-scale. Rock exhibits obvious nonlinear characteristics on the macroscopic scale, but from the microscopic point of view, the fracture properties of the local mesoscopic unit can be assumed to reflect elastic-brittle behavior. Based on this, it is appropriate to use the elastic-brittle constitutive relation to describe thermo-solid coupling deformation behavior at the mesoscopic level. Therefore, the above equations are still valid for this meso-scale numerical model.

Model configuration
The present study focuses on formation of surface cracks caused by cooling shrinkage. Due to the influence of natural environment such as weathering, rain erosion and temperature change, the mechanical properties of rock surface materials are generally different from interior materials. Therefore, we divide the model into two materials for analysis as shown in Fig. 3. The size of the model consisting of two layers is 150 mm × 150 mm × 30 mm. The upper layer is weak with a thickness h of 3 mm. The underlying material is substrate with a thickness of 27 mm. The model is divided into 5,400,000 hexahedral elements. To observe the development of cracks in the upper layer, the lower layer is made stronger than the upper layer. It is assumed that there is ideal contact between the two layers and no interface transition layer is inserted. The model is circumferentially adiabatic with a fixed temperature of 20 • C at the bottom and an initial temperature of 20 • C on the top. The cooling process is simulated by linear reduction of the surface temperature, and the temperature drops by 0.2 • C per step.

Digital image processing
From Fig. 1, it can be seen that the crack networks are intertwined, and the cracks and the rock blocks are the main components. For different sites, the fracture networks show significantly distinct geometric structure characteristics. The simulated result diagrams have to be processed for further quantitative analysis, as shown in Fig. 4. During digital image processing, the gray threshold value is determined to binarize the images using threshold segmentation technology. Then, the cracks and blocks can be represented by black and white pixels as shown in Fig. 4(b). In this way, the morphological data can be retrieved for analysis.
For comparing the morphological characteristics of the cracks and blocks, the measurement required are: (1) the surface crack rate, defined as the ratio of damaged area to total area. This reflects the cracking degree of rock; (2) the number of cracks; (3) the crack length; (4) the number of blocks and distribution of block area; (5) the fractal dimension of the crack network; (6) the acute angle between two cracks and (7) the number of polygon edges.

Analysis conditions
In order to study the influencing factors of polygonal cracks on the rock surface, the thickness of the upper layer, the rock heterogeneity and the boundary conditions, are varied. In this study, the total thickness of the model is maintained at 30 mm. The material parameters are shown in Table 1. The thickness h of the upper layer is set at 1 mm, 3 mm, 5 mm, 7 mm and 9 mm. In the simulation, the materials are assumed to be inhomogeneous, and the non-uniform distribution of the mechanical parameters obeys the Weibull distribution. The heterogeneity index m is used to characterize the degree of uniformity. The higher m is, the more homogeneous the rock is. In this study, the heterogeneity index m of the upper layer is set to be 2, 4, 6, 8 and 10. In addition, three kinds of boundary conditions are considered, i.e., biaxial constraint, uniaxial constraint and free boundary, as shown in Fig. 5. The formation, propagation and saturation of cracks was simulated under cooling.

Initiation and propagation of polygonal cracks
In this section, we describe the crack formation and propagation process of the surface layer as reproduced numerically under the biaxial constraints. The material parameters are listed in Table 1. The thickness of the layer h was set to be 3 mm and the model was circumferentially adiabatic. Meanwhile, a constant temperature of 20 • C was applied to the base. The surface temperature is represented by T, and initially a T of 20 • C was applied to the top. This was reduced by intervals of 0.2 • C per step as cooling proceeded.
From Fig. 6, it can be seen that the decrease of surface temperature leads to shrinkage effect of the surface layer. This shrinkage is constrained by the lower layer and boundary constraints, forming high tensile stress concentrations at the internal defects of the rock. As temperature T decreases, the accumulated tensile stresses continuously increase. Due to the inhomogeneity of the rock, the strength of each mesoscopic element varies. Low-strength elements damage first. At the initial stage of cooling, damaged elements are scattered randomly on the surface layer ( Fig. 6(b)). With the further decrease in T, more and more microcracks occur and begin to aggregate and connect, forming the identifiable cracks (Fig. 6(c)). At this time, even if temperature T drops only slightly, high stresses concentrate at the tips of cracks and promote   further crack expansion and formation of a macroscopic cracking network ( Fig. 6(d)). As shown in Fig. 6(e), many large polygonal areas form. However, as the temperature drops, new cracks generate inside the large polygons and gradually extend to the outer edges. This divides the large polygons into smaller ones. This insertion process of cracks is induced because the accumulated tensile stresses are still high enough to cause failure of new elements and a new round of crack propagation and coalescence, thus splitting the large polygons. It is worth noting that when the crack insertion process develops to a certain extent, even if the surface temperature is reduced again, it is difficult for new cracks to form and no new polygons are generated. This phenomenon is termed the crack saturation and has been clearly reproduced in Fig. 6(f). Fig. 7 shows the temperature field at the crack saturation stage, from which it can be seen that when there has been only a 39.8 • C drop in the surface temperature from 20 • C to − 19.8 • C, many surface cracks are produced. This emphasises that rock failure caused by thermal stresses should not be neglected in engineering situations. Note that in the FEMbased RFPA method, the small deformation assumption should be satisfied, and the initiation and propagation of cracks occurring in the small deformation state can be characterized by element damage. Therefore, the detaching, sliding and locking along the interface between blocks cannot be modelled. Meanwhile, the material is isotropic in current study.

Influence of the thickness of surface layer under biaxial constraint
To study the influence of the thickness of the surface layer on the formation of rock surface cracks, the thickness of the surface layer was set at h = 1 mm, 3 mm, 5 mm, 7 mm and 9 mm. The material parameters listed in Table 1 were used. Fig. 8 shows the simulated results of polygonal cracking of the rock surface with different thicknesses under the effect of cooling when the heterogeneity index m is 6.
It can be seen from Fig. 8 that, as the thickness of the surface layer increases, the average area of polygonal blocks increases, the total number of blocks decreases, and the shape of blocks becomes more regular. This results because the constraining effect of the lower base and boundaries is more significant, and the distribution of high tensile stresses is more dispersed than for a thinner surface layer under cooling shrinkage. Thus, the thinner the surface layer, the smaller the average area of an individual block, and the larger the number of blocks. In this way, combined with quantitative analysis, the thickness of a weak layer can potentially be calculated from the fracturing of rock surface.
It can be seen from Fig. 9(a) that as the thickness of the surface layer increases, the number of polygonal blocks formed by the cooling effect decreases. Clearly, when the thickness h increases from 1 mm to 5 mm, the number of blocks and cracks decreases. When the thickness h rises from 5 mm to 9 mm, these measures change only slightly. When the thickness of the surface layer grows to a certain degree, the number of blocks and cracks is insensitive to h. Meanwhile, the generation of cracks on surface layer is affected not only by thermal stress but also by constraints between the two layers. The binding ability between the lower base and the surface layer weakens as the thickness h increases.
The surface crack rate is defined by the ratio of the fractured area to the total area. The higher the surface crack rate, the greater the degree of block cracking. From Fig. 9(b), it is clear that as the thickness of the weak layer rises, the crack rate of the rock surface drops greatly during the first phase but tends to remain the same during the second phase. Therefore, thickness h has great influence on the fracturing patterns of polygonal cracks on a rock surface under the action of temperature stress when it changes within a certain range. The thicker the weak layer, the larger the average area of blocks is and the smaller the number of cracks.
The fractal dimension represents the complexity of the cracking network. The larger the fractal dimension, the more complex the network. As shown in Fig. 9(c), when the thickness of the surface layer is 1 mm, the fractal dimension is largest because the number of fractured blocks is largest. Hence, the crack structure is the most complex for this case. When the thickness h reaches 5 mm, the fractal dimension is smallest, indicating that the complexity degree of the cracking network is lowest. After that, the fractal dimension gradually rises at a small rate because of rock heterogeneity. Fig. 10(a) illustrates the distribution of polygonal block areas under different surface layer thicknesses. Most of the block areas lie in the range 0-22 × 10 3 pixels. As surface layer thickness increases, the proportion of large-area blocks gradually increases. Meanwhile, the number   of small-area blocks is obviously greater than large blocks when h = 1 mm, 3 mm and 5 mm. When h = 1 mm, 94% of the block areas are distributed in the range 0-2 × 10 3 pixels. Fig. 8(a) shows the distribution of small blocks generated when h = 1 mm. When h = 3 mm, 34%, 28% and 18% of the block areas are distributed between 0 and 2 × 10 3 pixels, 2-4 × 10 3 and 4-6 × 10 3 pixels, respectively. 80% of the block areas are distributed between 0 and 6 × 10 3 pixels, indicating that small blocks comprise the majority. Large block areas are mainly distributed in the range of 14-16 × 10 3 pixels.
When h = 5 mm, the distribution of block areas tends to be uniform, and the main range is 0-22 × 10 3 pixels. 47% and 40% are distributed between 0 and 6 × 10 3 pixels and 6-16 × 10 3 pixels, respectively. At this time, some large blocks begin to appear and account for 13% between 16 and 22 × 10 3 pixels. When h = 7 mm, 46%, 26%, 16% and 12% are distributed between 0 and 6 × 10 3 pixels, 6-16 × 10 3 pixels, 16-22 × 10 3 pixels and 22-34 × 10 3 pixels, respectively. When h = 9 mm, 30%, 40%, 14% and 5% are distributed in 0-6 × 10 3 pixels, 6-16 × 10 3 pixels, 16-22 × 10 3 pixels and 22-34 × 10 3 pixels and 38-40 × 10 3 pixels, respectively. Fig. 8(e) shows that the largest blocks appear in the model with the thickest weak layer. As the thickness of the weak layer increases, the number of small blocks decreases, and the number of large blocks increases. However, the number of small blocks still greatly exceeds that of the large blocks in terms of area distribution. Fig. 10(b) shows the distribution of all crack lengths for the models under different weak layer thicknesses. No matter how the thickness of the weak layer changes, the cracks with length 20-40 pixels comprise the majority. Above 40 pixels, the number of cracks gradually decreases as the interval value increases. Regardless of the change in layer thickness h, the number of short cracks generally exceeds that of long cracks. Over 60% of the crack lengths are concentrated in the range 10-70 pixels.
Because of the influence of polygonal area and shape, the length of a crack is limited and is related to the morphological characteristics of polygons. The numbers of edges of polygonal blocks with different weak layer thicknesses are shown in Fig. 10(c). With layer thickness increasing, the number of polygonal edges tends to be stable. Clearly, the number of quadrilateral, pentagonal and hexagonal blocks is obviously more than that of other blocks. This reflects what is directly observed in nature. From Fig. 10(d), the acute angle distribution between two cracks under different weak layer thicknesses can be analyzed. No matter how the thickness of the weak layer changes, acute angles in the range of 60 • -180 • comprise the majority. There are few acute angles between 0 • and 60 • . Additionally, the probability density of block area, crack length, polygon side and intersection angle can provide necessary data for generating discrete fractures and blocks in DFN analyses.
A new crack usually generates at pre-existing defects in rock.
Propagation starts from a damage point and develops along a path induced by the process of stress buildup, stress shadow and stress transfer. According to the minimum energy and shortest path principles, trifurcation appears, and the acute angle between cracks is generally about 90 • or 120 • , which has been simulated as shown in Fig. 11(a). Furthermore, the actual value is affected by the properties of rock materials. When the material is isotropic, once the fracture is initiated on the rock surface, it is common to trifurcate with 120 • intersecting angles. Simultaneously, when a major crack has already formed, the tensile stress perpendicular to the crack becomes dominant under loading. Therefore, if a new fracture intersects the existing crack, the intersection angle will be perpendicular, i.e., a 90 • intersecting crack will occur. From Figs. 6 and 11, we can see that our simulation has produced the progressive formation of the observed 90 • and 120 • intersecting cracks in field by characterizing the material heterogeneity and calculating the process of stress buildup, stress shadow and stress transfer. Fig. 12 shows the simulated polygonal cracks of the surface layer with different heterogeneity indices under cooling when the thickness h is 3 mm. It can be seen from Fig. 12 that the more uniform the rock is, the more regular the polygonal blocks are. When the heterogeneity index m = 2, there are mostly small, disordered cracks in the model. Although the number of cracks is large, the number of penetrating cracks is relatively small. With further growth of rock uniformity, the number of non-penetrating cracks reduces, as does the complexity of crack. Therefore, the heterogeneity of rock materials has a great influence on the shape, area, and distribution of polygonal blocks under the action of thermal stress.

Influence of heterogeneity of the surface layer under biaxial constraints
As shown in Fig. 13(a), as rock homogeneity increases, the number of polygonal blocks generated on the rock surface increases, and then gradually stabilizes after m = 6. When the heterogeneity index m = 2, although the number of cracks is not significantly different from other cases, the randomness of the failure site increases such that fine cracks cannot easily penetrate and irregular, large blocks form. With the increase in rock homogeneity failure becomes more regular, resulting in relatively regular blocks. Considering that when the heterogeneity index m = 2, the number of blocks is much smaller than for other cases, but the number of cracks is similar, it is clear that many cracks do not penetrate because of material heterogeneity. With the increase of rock homogeneity, the number of cracks tends to stabilize, and the overall trend is similar to the number of blocks.
It can be seen from Fig. 13(b) that, with the gradual increase of homogeneity, the surface crack rate of the rock increases but the numbers of cracks do not. This indicates that the number of penetrating cracks increases. The more uniform the rock is, the easier it is to form a penetrating crack. Fig. 13(b) demonstrates that the rock surface crack rate is an important indicator of the degree of uniformity of the rock. From Fig. 13(c), we can see that with the increase in homogeneity, the overall trend in fractal dimension is an increase. This is because the proliferation of penetrating cracks leads to more complex crack structure. However, the details are also influenced by the material heterogeneity.   Fig. 14(a) shows the distribution of polygon areas under different heterogeneity indices. Regardless of homogenization, more than 50% of the block areas is distributed in the range of 0-6 × 10 3 pixels. Except for the case of m = 2, the other block areas are mainly concentrated in the range 2-4 × 10 3 pixels. When m = 2, the area distribution of polygons is relatively dispersed, and more than 50% of block areas are concentrated in the range of 0-8 × 10 3 pixels. Although some large polygons exist, there are still more small blocks in general. When m = 4, 6, 8 and 10, a basically similar distribution law of block areas is found, and the block area is mainly distributed in the range 0-6 × 10 3 pixels. Due to the randomness of the numerical simulation, the heterogeneity index m = 6 varies slightly but the overall trend is consistent. During the process of crack formation and propagation, new cracks easily generate in largearea blocks, resulting in large blocks being divided into multiple small blocks. Thus, the number of small blocks generally exceeds that of large blocks.
As shown in Fig. 14(b), the number of cracks in the range 20-50 pixels is much higher than in other regions no matter how the homogeneity changes. Above 40 pixels, the number of cracks gradually decreases with the growth of the interval value. Regardless of the change of the rock homogeneity, the number of short cracks is significantly more than that oif long cracks because of the formation mechanism of polygonal blocks and the properties of rock material. Fig. 14(c) shows that regardless of rock homogeneity, the quadrilateral, pentagonal and hexagonal blocks are always the majority. Compared with thickness, the edge number of polygonal blocks is more sensitive to the heterogeneity of the rock material. When the rock material is more inhomogeneous, the polygonal shapes formed by fracturing are more complex. Ideally, the perfect hexagonal block is likely to form in a completely uniform material. Besides, Fig. 14(d) shows that no matter how the homogeneity changes, the number of crack angles in the range 60 • -180 • is larger than in other ranges. Only a few acute angles lie in the range of 0 • -60 • , and most are in the range 90 • -150 • . This agrees with the distribution rules of cracks for different thicknesses of the weak layer. In addition, the probability density of block area, crack length, polygon side and intersection angle can help to generate discrete fractures and blocks in DFN analyses.

Influence of different boundary constraints
In nature, the formation of polygonal cracks on a rock surface is also affected by different constraints. In this section, the origin and propagation of cracks were reproduced under different boundary conditions including uniaxial constraint and free boundary (circumferentially unconstrained). The material parameters listed in Table 1 were used. The thickness h of the weak layer was set to be 3 mm, and the model was circumferentially adiabatic with a constant temperature of 20 • C at the bottom and an initial temperature of 20 • C on the top. The surface cooling process was simulated by linearly reducing the surface temperature by 0.2 • C per step.

Uniaxial constraint
As shown in Fig. 5, the boundary condition is changed from biaxial constraint to uniaxial constraint along the Y axis. It can be seen from Fig. 15 that during the cooling process, the rock surface still experiences a strong shrinkage effect because of cooling. However, due to the uniaxial constraint, the tensile stress generated on the rock surface is mainly perpendicular to the Y direction. Thus, the cracks manly expand horizontally along the X axis. Similar to the biaxial constraint, with continuous cooling, the tensile stress accumulated on the rock surface grows as a result of thermal contraction. Because of the inhomogeneity of the rock, the distribution of element strengths is not uniform, and the cracks initially appear at the low-strength elements. Therefore, at the initial stage of cooling, surface microcracks are scattered randomly on the rock surface ( Fig. 15(b)). With further cooling progressively more microcracks occur, aggregate and penetrate, forming identifiable cracks (Fig. 15(c)). The formation of these cracks is a prerequisite for stress concentration. After that, even if cooling is not too great, high stress concentrations will form at the crack tips, promoting further expansion of the cracks. Most of the long cracks develop towards the left and right sides of the mode along the X axis, as shown in Fig. 15(d).
From Fig. 15(e), it can also be seen that the interval between initial cracks is large. However, as the temperature decreases, more small cracks occur between the pre-formed long cracks. Being restricted along the Y axis, tensile stress is mainly perpendicular to the Y axis, resulting in a series of parallel cracks along the X axis and the formation of a special cracking structure, as shown in Fig. 15(e). The cracks are not completely parallel due to the influence of material heterogeneity. When the crack insertion process achieves a certain extent, even if the surface temperature continues to reduce, new cracks are difficult to form, i.e., cracking is saturated, as shown in Fig. 15(f).

Free boundary
When there are no displacement constraints applied to the model, Fig. 15. Crack formation and propagation process on cooling rock surface under uniaxial constraint by the maximum principal stress contours (unit: MPa).
the rock surface is only affected by thermal stress, and the energy required for failure increases. As shown in Fig. 16(b), when the temperature drops below − 20 • C, the model surface begins to be damaged. Due to the influence of rock inhomogeneity, at the initial stage of cooling, the damaged areas are scattered randomly but relatively concentrated at the center of the model (Fig. 16(d)). During cooling, the concentration of high stresses at the tips leads the cracks to develop and connect. The development direction of the cracks was synchronously towards the X axis and the Y axis, forming several polygonal blocks at the center of the model, as shown in Fig. 16(e). When there are no displacement constraints applied to the model, the rock surface is only affected by thermal stress, and the energy required for failure increases. As shown in Fig. 16(b), when the temperature drops below − 20 • C, the model surface begins to be damaged. Due to the influence of rock inhomogeneity, at the initial stage of cooling, the damaged areas are scattered randomly but relatively concentrated at the center of the model (Fig. 16(d)). During cooling, the concentration of high stresses at the tips leads the cracks to develop and Fig. 16. Crack formation and propagation process on cooling rock surface under free boundary by the maximum principal stress contours (unit: MPa).
connect. The development direction of the cracks was synchronously towards the X axis and the Y axis, forming several polygonal blocks at the center of the model, as shown in Fig. 16(e). Similar to the constrained boundary conditions, the initial cracks form large polygonal areas. However, with cooling, new cracks generate in the polygons and gradually expand to the edges of the polygons, thus splitting the large polygons into smaller polygons, as shown in Fig. 16(g). But the difference between the biaxial-and uniaxial constraint is that when the cracks extend outward they only develop perpendicular to the model boundaries. Finally, the center part of the model surface presents polygonal cracks, and the areas near the upper and lower boundaries produce parallel cracks, as shown in Fig. 16(h). This phenomenon results from high tensile stress near the model boundaries being perpendicular to the boundaries. Meanwhile, cracking saturation also occurs as shown in Fig. 16(i).
These results explain the regularity and complexity of polygonal cracks on rock surfaces under cooling in nature to a certain degree. Meanwhile, this study provides a possible way for evaluating the upperlayer thickness, heterogeneity and historical stress environment with observed polygonal cracks in field. In terms of the shortcomings, note that the ratios of the physical and mechanical parameters between the surface layer and the substrate, such as elastic modulus and compressive strength, are, limited in the current study. However, such ratios may greatly influence the decay patterns and mechanism of rock surfaces (García-Rodríguez et al., 2015a, 2015b, which needs further study.
Additionally, the quantitative analysis of geometric features of polygonal cracks on rock surfaces with different thicknesses, heterogeneity indices and boundary constraints, including the probability density of block area, polygon side number, crack length and intersection angle, can help to determine the potential occurrence and shape of fractures for establishing DFNs in geo-hydrology. The in-depth mechanical calculation can be conducted accordingly.

Summary and conclusions
Using meso-damage mechanics, continuum mechanics and thermodynamics, we studied the formation and evolution of polygonal cracks on rock surfaces in nature under cooling. A heterogeneous model reflecting the characteristics of a surface weak layer was tested using the RFPA3D method. The morphological characteristics and mechanical mechanism of surface cracking was discussed. The main conclusions are as follows: (1) At the initial stage of cooling, damage points on the surface are scattered randomly. Then, damage begins to aggregate and connect, forming identifiable microcracks. With cooling these microcracks continue to expand and coalesce, forming macroscopic cracks. Meanwhile, the evolution mode of cracks changes under different boundary conditions. Under biaxial constraints, the first cracks generally form a large-area polygon. Due to cooling new cracks occur in the polygon area and gradually expand outwards. Thus, large polygons are gradually divided into small polygons. However, this insertion process of cracks cannot continue when the surface temperature is reduced more than a certain amount, i.e., crack saturation appears.
(2) The polygonal cracks or cracks parallel to the free direction form under biaxial-or uniaxial constraints. However, the rock surface presents polygonal cracks at the center and parallel cracks in the surrounding areas under unconstrained conditions. Under strong biaxial constraints, stress is basically isotropic and an intersecting angle of 90 • or 120 • is easy to form. This leads to the development of polygonal cracks. Under uniaxial constraints, the material cannot shrink freely along the constrained axis. Thus, high tensile stress results in a group of parallel cracks which are perpendicular to the constrained axis. For a free boundary, the energy required for failure accumulates at the center of the model, where polygonal cracks occur and expand outwards. Then, each crack develops along the direction perpendicular to the corresponding boundary because of the high tensile stress. (3) The thinner the weak layer, the smaller the average area of blocks. The largest blocks appear in the model with the thickest weak layer. As layer thickness increases, the number of small blocks decreases, and the number of large blocks increases. However, the number of small blocks is generally greater than large blocks because large blocks are more likely to produce inner tensile stress concentrations due to material inhomogeneity. They promote the generation of new cracks and divide large blocks into many small blocks. There are more short cracks than long cracks because of small polygonal blocks. (4) As rock uniformity increases, the number of non-penetrating cracks reduces, as does the complexity of crack shape. When the material becomes inhomogeneous, the randomness of damage increases, and irregular large-area blocks tend to form. When rock uniformity increases, the surface crack rate rises faster than the number of cracks, demonstrating that the number of penetrating cracks increases. Furthermore, with the increase of material homogeneity, fractal dimension also increases because more through-cracks lead to a more complex distribution of cracks on the rock surface. (5) The acute angles between two cracks mostly lie in the range of 90 • -150 • regardless of the change in layer thickness and material inhomogeneity. According to the minimum energy and shortest path principles, the acute angle is generally about 90 • or 120 • . However, the actual value is affected by the rock inhomogeneity. On the one hand, trifurcation tends to appear with 120 • intersections. On the other hand, if a major crack already exists, a developing crack will change its path to the direction perpendicular to the existing crack due to the high tensile stress parallel to the existing crack. Hence, a 90 • angle forms. Additionally, the edge number of polygon blocks is sensitive to rock heterogeneity, i.e., polygons tends to become more complex for a more inhomogeneous layer. Additionally, in terms of the computational requirements for practical DFN analyses, a rock with 0.001 m -10 m in length, width and height can be considered if the fractures are generated using the RFPA3D method.

Disclosure statement
The authors declare that they have no known competing financial interests or personal relationships that could appear to influence the work reported here.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The datasets generated and/or analyzed during the current study are available from the corresponding author upon reasonable request.