Stress interactions in cracked media during the closure of prestressed cracks

With the increasing pressure, a crack in a medium wi l l be gradually closed, which is affected by stress interactions. The closing process of parallel cracks under vertical stress is simulated here. The coplanar and stacked cracked models are constructed to analyze the influence of two types of stress interaction on the closing process. The spatial distribution of cracks, demonstrated by numerical experiments, has a significant impact on stress interactions and thus the process of crack closure. The mechanisms underlying the delay of crack closure caused by stress interactions are different for the two models. Furthermore, according to the stress dependence of crack microscopic parameters (crack porosity, aperture, and length of major axis), the process of crack closure can be divided into three stages: the linear deformation stage, the contact stage, and the closure stage. In the first stage, no contact is permitted inside the crack. The shielding effect directly leads to a closure lag, and thus a linear stress dependence of the microscopic crack parameter. In the second stage, the shielding effect determines the increasing rate for the stress dependence of microscopic crack parameters in regularly distributed cracked models. However, for the randomly distributed cracked model, local stress interactions result in the eccentric closure of cracks, and thus the crack closure lag. In the last stage, the crack is closed and stress interactions disappear both in the regularly and randomly distributed models.


Introduction
Due to the influence of tectonic movement and other factors, cracks are widely distributed underground in natural rocks (Cao et al. 2019, Wang & Lei 2021, Fu et al. 2022 ).The cracks wi l l affect the mechanical properties, permeability, and weathering resistance of the rock (Liao et al. 2021, Sotelo et al. 2021 ).In the process of flourishing oil and gas production, many unconventional resources exploration activities, such as dri l ling, oi l recovery, etc., have a significant impact on the stress distribution inside the rock, leading to deformation variations of the crack, and thus the elastic properties of the rock.Therefore, it is of great significance to analyze the effect of stress on crack to improve oi l wel l productivity, optimize development strategies, and protect reservoirs.
Generally, the stress dependence of the elastic modulus of the rock is controlled by the deformation of the soft crack, that is, cracks with small aspect ratio or contact areas between particles (Walsh & Grosenbaugh 1979, Deng et al. 2015 ).Current study in this field could be divided into three aspects: The first type is the elastic acoustoelastic theory of prestressed medium (Thurston & Brugger 1964, Sinha 1982 ), which has been developed in the field of metal materials.Based on continuum mechanics of finite deformation, this theory describes the acoustic-elastic effect caused by infinite deformation with small amplitude fluctuation, superimposed on the finite deformation of prestressed background (Pao 1984, Norris et al. 2007, Fu et al. 2020 ).The second type generally believes that soft cracks contribute to the nonlinearity of elastic properties (Gurevich et al. 2011, Collet et al. 2014, Glubokovskikh et al. 2016, Ren et al. 2021 ).It is assumed that the stress dependence of the contact area and pressure inside the crack follows the exponential law.Combined with the effective medium theory, the anisotropic parameters of the model can be estimated.The third type is the double-porosity medium model and its derivative models.The pores are divided into stiff pores with large aspect ratio and compliant cracks with smaller aspect ratio ( < 0.01).The ratio directly determines the nonlinear stress dependence of cracks.
A uniform stress distribution on the crack surface is generally assumed in most of these theories.However, previous studies have demonstrated that the stress distribution is actually nonuniform due to the heterogeneities of underground rocks, which brings a huge challenge to the applicability of these theories.
Studies have shown that when the distance between cracks is small, the stress interactions wi l l also significantly affect the spatial distribution of the stress on the crack surface.For example, the spatial distribution of cracks wi l l greatly affect the stress interactions (Grechka 2007 ).Based on this, the stress interactions are divided into the amplification and shielding effects, which indicate the areas with the larger and smaller stresses, respectively (Zhao et al. 2015 ).Generally speaking, the stress shielding area mainly occurs in the stacked cracked model, while the stress amplification area mainly occurs in the coplanar cracked model.Both stress interactions wi l l significantly change the spatial distribution of the stress inside the model, and thus further affect the macroscopic parameters of the model (Zhao et al. 2015 ).
For further i l lustrating the influence of two types of stress interaction on the closure of cracks, Wang et al. (2016 ) set up a 2D cracked model to systematically study the stress dependence of microstructure parameters during the crack closure.The results showed that when the crack density increases, the interactions between cracks wi l l significantly affect the spatial distribution of the stress on the crack surface (Cao et al. 2021 ) and reshape the process of crack closure, resulting in significant differences between the predicted closure stresses based on the classical model and the real ones.
However, most existing theories cannot distinguish the two stress interactions from the stress-dependence properties, leading to the loss of correlations between macroscopic properties and stress interactions.Furthermore, most previous studies (Cao et al. 2019 ) about the stress interactions assume the infinite deformation, which, however, are not suitable to describe the deformation during the crack closure.Therefore, the influence of two stress interactions on the crack closure has not yet been correctly clarified.
Motivated by this problem, the influence of the stress shielding and amplification on the crack closure process wi l l be numerically investigated.This paper is organized as follows.Section 2 introduces the chart flow for deriving the microscopic crack parameters and compressive modulus in each vertical stress step.Then, two typical models (the stacked cracked and coplanar cracked models) are used to study the characteristics of the stress interactions during the crack closure.Finally, effects of stress interactions as well as the stress dependence of microscopic cracked parameters are studied in the regularly and randomly distributed cracked models, while considering the spatial distribution of cracks.

Methodology description
We constructed a 2D rock sample containing one elliptical crack to i l lustrate the numerical scheme (Fig. 1 ).
The modeling sample consists of: a linearly elastic solid matrix and an elliptical crack (Fig. 1 ).The stress-strain relations of the rock sample satisfy the following equation: (1) In this equation,  and  are the stress and strain tensors, respectively.C is the elastic stiffness tensor.For an 2D isotropic sample, the vertical component of the stress  yy is defined as where  and  are the lame coefficients of the matrix,  yy is the vertical strain.Therefore, we could use the  yy and strain  yy to calculate the compressive modulus (P-wave modulus, H ), which is defined as For a case with the boundary condition set (  xx = 0), the Pwave modulus ( H ) can be determined from the stress-strain ratio based on Equations ( 2) and ( 3): The load F is applied on the top of the model: where n is the unit boundary normal, F is the force per unit area on the top boundary, p is the pressure on the top boundary, and y 0 is the position where the loading is applied (top boundary): where u = ( u x , u y ) is the displacement vector.n r u = 0, means that roller boundary is used at the lateral boundaries, that is, the normal displacement at the boundary is zero.We assume the stress-free boundary and fixed boundary conditions for surfaces of the crack and bottom of the model, respectively.
According to the aforementioned formula, the displacement field u( x, y ) of the crack inner wall can be calculated first.Then, based on the displacement, the corresponding crack deformation can be obtained.

Calculation of the crack deformation
Under each load step, the deformation for the inner boundary of crack is computed to determine closing process of the crack.If the y -coordinates of points on the inner wall of the crack are identical, the crack wi l l be closed.However, if the crack is not closed, the geometric cracked parameters are updated according to the displacement field.Then we continue to increase the compressive load under Δ i .This process wi l l cycle until complete closure of all cracks.
The length ( d ) of major axis for the elliptical crack under compression was previously assumed to be stressindependent (Wang et al. 2016, Ren et al. 2020 ).However, the stress dependence of d inside the crack has been verified by the compressive test (Sarout et al. 2017 ).Therefore, besides the micro-crack parameters, such as aperture and porosity, the stress dependence of d is also systematically investigated to reveal the crack closure process.
Parameter 1 (crack aperture): The y -coordinates of each point in the upper and lower boundary of the crack, are updated according to the displacement of crack boundaries in each compression step [Equation ( 7)].Then, the aperture of the crack [ Aperture in Equation ( 8)] is computed as the difference of y -coordinates between the upper and lower surfaces [Equation ( 8)] where the superscript i represents the index of the loading step, x and y are the x -and y -coordinates of the inner wall of the crack, respectively, and v is the displacement of the points in the crack wall.The subscript of lower and upper represent points at the lower and upper boundaries of the crack, respectively.
Parameter 2 (crack porosity): The porosity is the area of the deformed crack at each loading step.Wang et al. (2016 ) adopt a series of rectangles to discrete the ellipse crack, which, however, is expensive in computation.To simplify the solution formula of area, we adopt the Gauss's theorem to convert the surface integral to the boundary integral, and thus calculate the area of the crack as where the superscript i represents the index of the loading step, x is the x -coordinate of the inner wall of the crack, and n x is the x component of the outward normal vector of the boundary, °C is the contour profile of the crack.Based on the crack area, Equation ( 10) is adopted to obtain the crack porosity of the model: where A is the model area.
Parameter 3 (averaged length ( d ) for major axis of crack): Based on the data from steps 1 and 2, we can calculate the coordinate of each point of all cracks after loading.
The two points [ x , y lower ( x )] and [ x , y upper ( x )] in the two surfaces of the crack wi l l be in contact if the y -coordinates are identical.Furthermore, two adjacent point contacts indicate the parts between them are in contact.Otherwise, the part is open.The coordinates of the crack surface and the major axis ( d at the bottom of Fig. 2 ) of the crack continue to be updated 157 Downloaded from https://academic.oup.com/jge/article/21/1/155/7499738 by guest on 26 September 2024 Considering stress interactions, each crack in the sample deforms differently.Therefore, the averaged crack parameter is highly appreciated to represent overall characteristics of all cracks for each loading step.Then, the averaged crack aperture, porosity, and major axis ( d ) can be calculated by Parameter 4 (static P-wave modulus): For a linearly elastic material, the vertical compression modulus H with the boundary condition in Equation ( 6), is defined as Here Δ i yy and Δ i yy represent changes in the vertical (axial) stress and strain variations in each loading step, respectively.The two parameters are defined as where u i y (x, y 0 ) is the vertical displacement of the top boundary in each loading step i , L i is the vertical length of the model in each compressive stress, and p i (x, y 0 ) is the pressure on the top boundary.
It should be noted that all the y-coordinates in Equations ( 7) and ( 8) correspond to the points at the crack walls.However, the strain in Equation ( 14) relates to the points at the top boundary.Therefore, the y -coordinates in Equations ( 7) and ( 8) are irrelevant to the strain in Equation ( 14).

The dynamic modeling scheme
The dynamic process about the numerical simulation is shown in Fig. 3 .The crack deformation is simulated under each load step  i .The crack structure parameters, such as aperture, porosity, and length of major axis for the crack can be updated based on Equation ( 11 ).As stress reaches the maximum value of the loading (  max ), our dynamic simulation stops, or the stress continues to increase with stress step Δ.

Model description
Assuming that the rock matrix is a homogeneous linear elastic body, a coplanar model (Fig. 4 a and b) and a stacked model (Fig. 4 c and d) are set up to quantitatively investigate the influence of stress interactions.The amplification and shielding effects dominate the coplanar cracked model and stacked cracked model, respectively (Cao et al. 2021 ).The distance between the cracks determines the degree of stress interactions.Therefore, the influence of stress interactions on the crack closure can be estimated by controlling the distance between cracks.
The model with a side length of 0.2 m, contains three elliptical cracks.The cracks have the lengths of 0.018 and 0.0002 m for major and minor semi-axes, respectively.Meanwhile, the model is assumed to be elastic during the deformation process.The elastic parameters for the matrix are given in Table 1 .

Mesh analysis
The zoom we selected, is a part of Model d (Fig. 4 d) containing three cracks inside.And the unstructured triangular mesh to discrete the zoom (Fig. 5 ), is used in such a way that the crack boundary coincides with the element boundaries.The largest and smallest mesh interval are about 0.004 m concentrated at the background area and 0.000013 m focusing at the tip of the crack, respectively.

Influences of stress interactions on crack closure
For better i l lustration of the influence of stress interactions upon the microscopic crack parameters during closure, a series of compression simulations are conducted: 5.7 and 11.7 MPa for coplanar models a and b (Fig. 6  Table 1.Elastic parameters for the matrix (Wang et al. 2016 ).

Elastic parameter Value
Density,  (g/cm 3 ) 2.02 P-wave velocity, V p (m s − 1 ) 2118.9 S-wave velocity, V s (m s − 1 ) 1254.7 Bulk modulus, K (GPa For the shielding effect on crack closure, comparison between two cases (Fig. 6 c and d  The numerical mesh to discrete the model is displayed.The selected numerical domain is a rectangular zone with three cracks inside (bold ellipse).For the mesh surrounding the crack, the spatial resolution of the mesh can vary strongly, ranging from 0.000013 to 0.004 m.
x -axis in Fig. 7 c).It is because the shielding effect greatly decreases the stress between cracks, leading to the lag in the crack closure.

Influences of stress interactions on static compressive modulus
As for the static compressive modulus (Fig. 8 ), it can be observed that with the increasing pressure, a slightly decreasing trend is dominated in the first stage.This is because compared with the background medium, the cracks are closed easier, leading to greater strain and thus smaller static compressive modulus of the model.In the second stage, more cracks begin to contact, resulting in a rapid increment in the static compressive modulus (Argatov 2021 ).
In the first stage, affected by the stress interactions, the order of modulus in the four models is Model b < Model a < Model c < Model d.It can be explained by that stress amplification makes the model "soft" while stress shielding makes the model "hard" (Zhao et al. 2015 ).The amplification effect dominates the coplanar model while the shielding effect dominates the stacked model.Therefore, Model b is the softest due to stress amplification, thus the cracks inside wi l l begin to be closed in advance.However, the transition pressure points from stress-dependence modulus to stressindependent modulus, are 20 MPa for both Model b and Model a (Fig. 8 ), indicating the stress amplification has little effect on the closure stress.The shielding effect of Model d is stronger that of Model c, resulting in the lag of crack closure.Thus, it is expected that the static compressive modulus (Model d) continues to increase beyond the stress of 60 MPa (pink line in Fig. 8 ).
Besides, an analytical solution (Guo et al. 2019 ) without the consideration of stress interactions, is also applied to the models presented in Fig. 4 to validate the accuracy of the numerical simulation.According to the characteristics of stress interactions, the influence of stress interactions is negligible for Model a in four models (Fig. 4 ).Our result shows that for Model a, the compressive modulus is 7.55 GPa by numerical simulation, and 7.76 GPa by analytical solution (The star * in Fig. 8 ).The difference between the two results is < 3%, verifying the small effect of stress interactions, so this is the accuracy of our numerical solution.

Model setup
The stress interactions are easily affected by the spatial distribution of the cracks.As a result, two types of cracked model with different spatial distribution are introduced to study the stress interactions caused by the spatial distribution of the cracks inside the models.Each square model has a side of 0.2 m and 20 ellipsoidal cracks.Each crack with the same aspect ratio, is organized in a different pattern.The dashed ellipsoid in Fig. 9 a represents the aspect ratio ( Ξ) of spatial distribution.In fact, the concept of the aspect ratio of spatial distribution represents the conditional probability of finding another inclusion given the position of an inclusion.A decrease in Ξ indicates the crack ends are approaching toward each other, suggesting an increase in amplification effect (Zhao et al. 2015 ).Conversely, increase in Ξ corresponds to increase in shielding effect.Therefore, the shielding effect is expected to be the greatest for Model 3 (Fig. 9 c), but the weakest for Model 1 (Fig.

Influences of stress interactions on crack porosity and aperture
As shown in Figs. 10 a and 10 b, with the increasing pressure, cracks in Models 1-3 wi l l be closed, and thus the porosity and averaged aperture of cracks wi l l gradual ly decrease.The process can be classified into three stages.The first is a linear deformation stage, where the porosity and averaged aperture of cracks are inversely proportional to the pressure.
In the second stage, the porosity and averaged aperture of cracks decrease at a smaller rate.In the third stage, the stressindependent parameters during the pressure range indicate the cracks are closed completely.As shown in Fig. 10 , all three stages can be observed in Model 1 and Model 4 (randomly distributed cracked model), while only the first and second stages can be observed in Model 2. Within the stress range of

Effect of stress interactions on length for major axis of the crack
As shown in Fig. 11 , the stress dependence of the averaged length ( d ) for major axis of the crack differs slightly from that of the crack porosity.All three stages are observed in Model 1.At the first stress stage within vertical stress range (0-10 MPa), the averaged length ( d ) for major axis of the crack remains unchanged, indicating the cracks do not contact.At the second stage within vertical stress range (10-20 MPa), the averaged length for major axis decreases rapidly with the increasing vertical stress.It can be explained that the contact area increases rapidly.At the last stage within the vertical stress range (20-60 MPa), the averaged length ( d ) for major axis continues to decrease at a negligible rate.Both the first and the second stages are observed in Model 2. However, only the first stage is observed in Model 3, indicating the 163 contact area is almost negligible.From Model 1 to Model 3, the stress shielding effect becomes stronger.Therefore, it can be inferred that the stress shielding effect delays contact and closure of cracks.
Comparing Figs 10 and 11 , it can be deduced that when the porosity and the aperture of cracks approach to zero, the averaged length for major axis of the crack is not zero, indicating that the crack has not completely been closed yet.

Effect of stress interactions on stress distribution on the crack
The stress distribution on the crack surface determines the geometric deformation of the crack on the microscopic scale, which wi l l be investigated further in the following.
It should be noted that since free boundary condition is set on the crack surface, the stress on the crack surface wi l l be zero.However, with the increasing vertical stress, the upper and lower boundaries begin to contact each other, leading to non-zero surface stresses, and the stress distribution density is not zero anymore.Therefore, the stress distribution density can be used to represent the contact var iations dur ing the closure process of cracks.
We believe that as the pressure approaches the closure stress, all cracks will be closed gradually.Furthermore, as the pressure grows continual ly, al l the microscopic crack parameters wi l l be stress-independent after the closure of the cracks.Therefore, if the parameter difference between the adjacent stresses is less than 1%, the stress is set as the closure stress.According to Fig. 10 , the closure stresses for Model 1, Model 2, and Model 4 are ∼20, 60, and 30 MPa, respectively.For Model 3, at the compression of 60 MPa, the microscopic parameters (porosity or aperture) have not reached the zero yet, meaning the closure stress is much more than 60 MPa.Therefore, deriving such an accurate closure stress wi l l result in expensive computation.For the simplicity, we set 60 MPa as the closure stress for Model 3. The stress distributions for  For Model 1, when the vertical stress approaches to 20 MPa, the stresses on crack surfaces are almost the same for all cracks (Fig. 12 a), indicating that all cracks wi l l be closed simultaneously.Moreover, the similar stress distribution on individual crack surface suggests that the negligible effect of stress interactions.
For Model 2, the stress shielding effect leads to smaller stresses in the zone between the adjacent cracks in the lower part of the model.Therefore, the cracks on the top are completely closed while other cracks are partially closing.
Due to the boundary conditions (Yurikov et al. 2017 ), a concentration of stress wi l l be induced inside the sample, which, similarly, is observed in the middle of the top cracks in Model 3 (black solid circle in Fig. 12 c).It should be noted that the abnormal stress distribution only exists in the middle of the top crack, leading to a negligible effect of the boundary conditions.Meanwhile, the stress shielding effect keeps most of the cracks open.It is consistent with the result in Fig. 11 that the averaged length ( d ) for a major axis of the crack remains unchanged.
It can also be observed from Fig. 12 d that the stress distribution on the crack surface is heterogeneous for the randomly distributed cracked model.On the contrary, the stress distribution on the crack surface is generally symmetrical for the regularly distributed cracked model.The stress interactions are believed to be the main factors inducing such a difference.

Influence of stress interactions upon the static compressive modulus
Generally, at the first stage, the static compressive modulus decreases weakly due to the crack closure without contact inside the crack.However, in the second stage, the modulus increases due to the crack contact.Finally, when the cracks are completely closed, the static compressive modulus of the model is equal to the background value.
For the regularly distributed cracked models ranging from Model 1 to Model 3, the shielding effect grows gradually, leading to the greatest static compressive modulus for Model 3 and the smallest for Model 1.In addition, with the increasing vertical stress, the shielding effect wi l l lead to a lag in the crack closure and thus a slower increase in the static compressive modulus.Therefore, the stressdependence range is the widest for Model 3 and the narrowest for Model 1.
Owing to the cancellation of the shielding and amplification effects, the static compressive modulus of Model 4 (randomly distributed cracked model) is similar to that of Model 1 at the first stage.However, the heterogeneous stress distribution of Model 4 directly leads to eccentric closure of cracks in the second stage (Fig. 13 ), thus a slower increase in the static compressive modulus compared with that of Model 1.

Discussion
The closure stress for the crack means that at this stress, both walls of the crack make contact with each other.The crack closure stresses for averaged crack porosity, av- eraged crack aperture, and macroscopic static compressive modulus are almost identical.For the elliptic crack with small 2D aspect ratio, the closure stress ( p closure ) can be calculated by (Mavko et al . 2009 ) where E m , v m , and  are the Young's modulus, Poisson ratio of the background matrix, and aspect ratio of the individual crack, respectively.According to the definition of closure stress, all the cracks wi l l be closed at the same stress if they have the same aspect ratio.According to Equation ( 15), the crack closure stress for the Models 1-4 in Fig. 6 is ∼22.9MPa.The simulated closure stress for Model 1 without stress interactions is ∼20 MPa (Fig. 10 ), verifying the reliability of our numerical simulation.
For randomly distributed cracked sample, the stress dependence of crack porosity of Model 4 (pink line in Fig. 10 a) is like that of Model 1 in the first stage while significantly different in the second stage.Such a difference indicates that the influence of stress interactions is different in the two stages.In the first stage, for Model 1, the stress distributions on the surface of all cracks are similar to each other (Fig. 10 a), suggesting that stress interactions have little effect on the deformation of cracks.Meanwhile, in the randomly distributed cracked sample, the stress interactions are also negligible due to cancellation of stress interactions (Zhao et al. 2015 ), making the stress dependence of physical parameters in Model 1 consistent with that of the Model 4 in the first stage.However, in the second stage, the symmetrical and eccentric closure (Fig. 12 ) patterns dominate for Model 1 and 4, respectively.The eccentric closure pattern wi l l produce a lag in the crack closure, and thus a wider range for stress dependence of microscopic parameters, such as the static compressive modulus.

Conclusion
In our study, based on the linear elastic hypothesis of the crack, vertical compression numerical tests are conducted.The stress dependence of the micro-geometric parameters of the crack and the macro-static compressive modulus are investigated during the crack closure process.Moreover, the influence of stress interactions on the closure process is studied.
Our results show that stress interactions (shielding and amplification) wi l l significantly change the crack closure process: amplification results in the eccentric closure for cracks leading to the closure lag, but has little effect on the closure stress (Fig. 8 ).On the other hand, shielding results in stress shadow areas between cracks.The stress shadow reduces the stress on the upper and lower surfaces of cracks thus significantly delays the closing process of cracks.
Based on the stress dependence of microscopic crack parameters, such as the aperture and the averaged major axis of the crack, the closure process can be divided into the linear deformation stage, the contact stage, and the closure stage, respectively.The stress interactions wi l l significantly affect the closing behavior of cracks.
In the first stage, the averaged crack aperture and crack porosity change linearly with pressure, while the major axis of the crack remains unchanged during this process, indicating that the upper and lower interfaces of the crack do not contact, and the static compressive modulus approximately remains unchanged.In the second stage, the averaged crack aperture and crack porosity change nonlinearly with pressure.The averaged length for major axis of the crack begins to decrease (Fig. 11 ), suggesting the increasing contact inside the crack, and thus a sharply increasing static compressive modulus.In the third stage, the crack is closed, and the averaged aperture, porosity, and major axis of the crack are approaching zero.Even when the static compressive modulus of the model approaches that of the background medium, the crack is not completely closed.
In the three stages, the crack interactions wi l l greatly reshape the spatial distribution of the stress on the crack surface, as well as the closure process.During the first and second stages, the interactions, especially the shielding effect, result in a smaller value of stress zone between the adjacent crack, which directly delays the closure of the crack.In the third stage, the stress interactions disappear.For the randomly distributed model, in the first stage, the two stress interactions can roughly cancel each other out.In the second stage, due to the effect of stress interactions in local scale, the eccentric closure of the cracks, directly leads to the crack closure lag.In the last stage, the crack is closed and the interactions disappear.

Figure 1 .
Figure 1.The 2D model of the cracked sample under compressive stress F .

Figure 2 .
Figure 2. Flowchart for calculating the dynamic averaged length ( d ) for a major axis during the crack closing process with increasing compression.

Figure 3 .
Figure 3. Flowchart for dynamic modeling with increasing compression.

Figure 4 .
Figure 4. Three cracked models.(a) and (b) Coplanar cracked models with different horizontal distances.(c) and (d) The stacked cracked models with different vertical distances.

M
= 2.02 × 2118.9 2 = 9.06 GPa For the amplification effect on crack closure, comparison between two cases (Fig. 6 a and b) suggests that the amplification intensifies the eccentric closure for the left/right crack, and accelerates the closure of the intermediate crack.For example, the profile of left crack in the model at the compres-sion of 11.7 MPa (solid red line in Fig. 7 a) shows no degree of eccentric trend, only a slightly eccentric trend for compression of 17.7 MPa (solid black line in Fig. 7 a).However, for the same crack in Fig. 7 b, the closure trend has already exhibited an eccentric tend at the compression of 11.7 MPa (solid red line in Fig. 7 b), and is more pronounced at 17.7 MPa (solid black line in Fig. 7 b).On the other side, the intermediate crack is sti l l open in Fig. 7 a (dashed black line in Fig. 7 a), but completely closed in Fig. 7 b (the dashed black line has coincided with the x -axis in Fig. 7 b).
) shows that the shielding effect delays the closure for all cracks, especially the intermediate one.For example, under the pressure of 53.7 MPa, the intermediate crack of the model in Fig. 7 d is sti l l open (black dashed line in Fig. 7 d), while the same crack in Fig. 7 c is completely closed (black dashed line has been coincided by the

Figure 5 .
Figure 5.The numerical mesh to discrete the model is displayed.The selected numerical domain is a rectangular zone with three cracks inside (bold ellipse).For the mesh surrounding the crack, the spatial resolution of the mesh can vary strongly, ranging from 0.000013 to 0.004 m.

Figure 6 .
Figure 6.Crack deformation for four models (Fig. 3 ) under various compression.For better visualization, horizontal and vertical scales are different in 6 (a)-(d), respectively, meanwhile all apertures are multiplied by 10.The solid line and the dashed line represent the contours of the cracks under the compression of 17.7 and 53.7 MPa, respectively.
9 a) as the normal stress is applied on the top of model.The finite element simulation is performed to get the crack closure of the randomly distributed cracked model.For the randomly distributed models (Fig. 9 d), the numerical models with a given crack density run 20 times by changing locations of the individual crack center.The corresponding results are collected for a further study.

Figure 7 .
Figure 7.The aperture distribution for four models, where (a)-(d) correspond to the four models (a)-(d) in fig. 4 , respectively.In the legend, LF, MF, and TF represent the left, middle, and top cracks, respectively.x is the major axis of the crack.The stress values in the legend correspond to different vertical stresses.

Figure 8 .
Figure 8.The stress dependence of the static compressive modulus.Models a-d correspond to the four models in Fig. 4 , respectively.The black label * corresponds the modulus without the compression, computed through the analytical solution.

Figure 9 .
Figure 9. 2D synthetic samples/blocks used in the analysis with 20 aligned cracks.The aspect ratio ( Ξ) of the crack distribution illustrated in (a).The spatial aspect ratio Ξ is the ratio between l and n ( Ξ = l/n ).Generally, the spatial aspect ratio increases from Model 1 to 3, corresponding to 1.05, 4.86, and 22.22, respectively.For Model 4, all the cracks are distributed randomly.

Figure 10 .
Figure 10.Stress dependence of crack porosity (a) and crack aperture (b).Models 1-4 correspond to Fig. 9 a-d, respectively.The centers of the bars for Model 4 are placed at the mean values of the numerically computed parameters (crack porosity and averaged aperture); the sizes of bars are their standard deviations.

Figure 11 .
Figure 11.Stress dependence of the averaged length ( d ) for major axis of cracks with the initial value of 0.036 m.Models 1-4 correspond to Fig. 9 ad, respectively.The centers of the bars in Model 4 are placed at the mean values of the numerically computed diameter; the sizes of bars are their standard deviations.

Figure 12 .
Figure 12.The stress distribution for the models (Fig. 9 ), where the length of the arrow represents the stress value.The vertical pressures are 20, 60, 60, and 30 MPa for the four figures, respectively.The arrow scale in black frame of the four subfigures, indicates the stress value of 100 MPa.The zoom-in circled by the black solid line in c is the area affected by the boundary condition.

Figure 13 .
Figure 13.The stress dependence of the static compressive modulus for Models 1-4 (the random model) in Fig. 9 .The centers of the bars in Model 4 are placed at the mean values of the numerically computed static compressive modulus; the sizes of bars are their standard deviations.