Numerical simulation of deformation memory effect of rock materials in low‐stress condition using discrete element method

The slip behavior of the pre‐existing sliding planes (eg, cracks or inter‐grain boundaries) in rocks is a natural candidate for a mechanism of the deformation memory effect. Based on this mechanism, the deformation rate analysis (DRA) method can be used to measure in situ stress in low‐stress region. However, the traditional discrete element method (DEM) is unable to model the deformation memory effect of rock when the loading stress is lower than the crack initiation stress due to the insufficient consideration for sliding planes. The simulated stress‐strain curves obtained from previous DEM are linear in initial loading stage until cracks start to grow. In the present work, the properties of planes (ie, cohesion, frictional resistance, and relocking mechanism) are introduced into DEM contact model to simulate the deformation memory effect of rock in low‐stress region, and its feasibility was confirmed by comparing it with experimental results. It is found that there exists a suitable in situ stress range for DRA method when using it in low‐stress region, among which the inflection is precise. Four factors related to the sliding planes are discussed with parametric study to study their influence on the deformation memory effect. The results show that the cohesion and dip angle are highly related to the threshold value, while the friction coefficient and fraction of slide planes have little influence on the deformation memory effect in rocks.

Among the existing measurement methods, 2 the DRA method proposed by Yamamoto et al 9 based on the rock deformation memory effect has attracted many practitioners. The DRA method is based on the compression tests in the laboratory on the in situ rock cores obtained from projects. In the past two decades, this method has been widely used in the field of rock engineering, such as reservoir construction 10 and tunneling engineering. [11][12][13] For underground resources production, since 2002, the DRA method has been used to measure the in situ stress in more than 40 areas of various types of mining projects in different countries such as Australia, Finland, the United States, and Canada. [14][15][16][17][18][19][20] The DRA method has been proved to be economical, efficient, and its accuracy can be improved by statistical analysis.
In brief, the DRA method is to apply several compression cyclic loads on the tested rock cores to determine the in situ stress. The cyclic loads consist of two parts, as shown in Figure 1A: (a) One is the preloading cycle which is used to produce "artificial stress" inside the rock sample to simulate the in situ stress and its peak loading stress is p ; (b) The other part is the measurement process containing two identical loading-unloading cycles (i cycle and j cycle), in which the applied maximum stress is m . After tests, the following equation is used to determine the in situ stress: where Δ i,j ( ) is the strain difference between two loading-unloading cycles, j ( ) and i ( ) represent the axial strain in the jth and ith stage, and is the corresponding axial stress. The typical DRA curves in axial and lateral directions are shown in Figure 1B. The DRA curve in the axial direction is concave in shape, and in the lateral direction is convex in shape. 18 The stress value in the x-axis of the inflection points corresponds to the maximum stress in the preloading stage, in other words, the in situ stress.
There are numerous researches on deformation memory effect and DRA method since it was proposed. Some researchers 19,[21][22][23] concluded that the mechanism of rock deformation memory effect was the same as the Kaiser effect 24 -the propagation of pre-existing cracks and generation of new cracks in compression. When the loading stress exceeds the crack initiation stress ci , the cracks will initiate and propagate. 25 The cracking behavior in a rock sample leads to nonlinear strains, which include both reversible and irreversible strain. The reversible ones consist of frictional sliding over the sliding plane, isolated tensile cracks opening and closure, and a change in the density of the tensile cracks. All this kind of nonlinear behavior in strain was considered to be reversible during the cyclic loading in the DRA test. Thus, the reversible strain can be canceled by Equation (1). The irreversible stain component which could not be canceled results from the growth of pre-existing cracks and the generation of new cracks when the applied stress exceeded the peak stress previously applied. Therefore, theoretically, a strain difference can be obtained when using Equation (1), and an inflection point will appear in DRA curve at the previous peak stress point. In this type of theory, as the deformation memory effect is due to the crack initiation and propagation, we referred this mechanism as the "cracking mechanism" (Figure 2A).
However, there are many experimental evidences indicating that the in situ stress, which is below the crack initiation threshold, can also be measured by DRA method. A limited summarize of strain difference-stress curves by experiment for various rocks is shown in Figure 3. Yamamoto et al 9 reported that for granodiorite core samples, the DRA method could detect the in situ stress state of approximately 1-6 MPa ( Figure 3A). Seto et al 19 and Hunt et al 26 showed that the DRA inflection existed when the prestresses were less than 15% of the UCS ( Figure 3B,C). Attar et al 11 found DRA inflections both in the brittle rock (Tuff) and soft rock (Sandstone) in the crack closure region ( Figure 3D). Chan 27 showed that the DRA method could determine the prestress that was less than 20% of UCS ( Figure 3E). Hsieh 28 conducted DRA test on sandstone and found that there is a kink in the DRA curve at approximately preloaded stress value, which is about 10% of the UCS ( Figure 3F). In addition, according to Yamshchikov et al, 23 the deformation memory effect existed in both elastic deformation and plastic deformation stages. Thus, further investigations are needed to explain the deformation memory effect when in situ stress is low.
Some researchers focused on the sliding occurred at grain boundaries and weak planes, and they thought that this was a natural candidate for the deformation memory effect in low-stress region. Stevens et al 29 concluded that the sliding crack model could reproduce the memory of the previous peak stress. Wang et al 30 proposed a "frictional sliding mechanism" for layered rocks attempting to explain the deformation memory effect in low-stress condition. It is well known that rocks are typical materials which contain a number of pre-existing cracks or interlayer boundaries and they are three-dimensional in nature, [31][32][33][34] and they are served as "sliding planes" in compression. Under loading, these planes will slide when the shear stresses exceed their sliding criterion sl . Subsequently, they would stop sliding when the shear stresses are smaller than the sliding criterion, and some residual elastic energy will store inside this slide system, which leads to the increase of sliding criterion. Thus, in the following loading stage, only when the shear stress exceeds the new sliding criterion, the sliding would occur again and lead to the change of slope of the stress-strain curve, and thus the rock memory phenomenon occurs ( Figure 2B). Hsieh 28 developed this model and made it not only suitable for layered rock but also other types of rock, which contain randomly distributed small pre-existing cracks or grain boundaries.
However, in order to reduce the number of parameters, there was an oversimplification in the proposed "friction sliding" model 28,30 that the influence of the compressive stress normal to the sliding planes was neglected. In addition, the used simulation tool could not highly replicate the rock property, such as the discontinuity and discrete properties. 35 In the present work, a commercial software Particle Flow Code (PFC) 36 based discrete element method (DEM) 37,38 was used which could highly simulate the features of rock materials. Some researchers modeled the DRA method based on PFC. 26,39,40 Hunt et al 26 concluded that if the loading stress was below the crack initiation stress the deformation memory effect would not be observed because there was no crack generated, which may be attributed to the fact that the traditional "crack model" was not suitable for this condition. Thus, the aim of this present work is to F I G U R E 2 Two mechanisms for deformation memory effect: (A) Crack mechanism; (B) Frictional sliding mechanism F I G U R E 3 Stress-strain difference curves for various rocks of experiments summarized from references 11,23,[26][27][28]51 3030 | TANG eT Al.
introduce a more realistic "friction sliding" model considering cohesion, frictional resistance, and relocking mechanism properties to simulate the deformation memory effect in low-stress region. The long-term aim is to numerically examine sample behavior with different rock types and properties under a variety of initial stress conditions based on the established model.
The structure of this paper is as follows: In Section 2, a methodology is introduced to model the deformation memory effect in low-stress region based on PFC3D; Section 3 provides a discussion of the results obtained and proves the feasibility of the model; In Section 4, the effect of main microparameters on deformation memory effect is investigated systematically, and these parameters are related to the slide planes, including cohesion, friction coefficient, dip angle, and fraction. Some suggestions for DRA method when facing low-stress condition are proposed. Finally, some conclusions are drawn in Section 5.

| Particle Flow Code
PFC is based on the principle of the discrete element method and can be used to study the mechanical behavior of rocks and rock-like materials (eg, granite and concrete). 41,42 In PFC, the material is modeled as collections of particles (rigid circular disks in 2D and sphere in 3D). Each particle is in contact with the neighboring particles via the contact element. As PFC2D assumes planar conditions, it cannot provide an accurate mechanical behavior of a cylindrical rock sample. 26 Thus, we adopted PFC3D in the present work.

| Modeling the frictional sliding mechanism
In previous DEM numerical simulations on DRA method, 26,39,40 when the applied stress is smaller than the crack initiation stress, there is no crack generation and propagation, so there was no "memory" occurred. The reason is that these studies insufficiently considered the characters of rock materials as discussed above.
As can be found in Wang et al's research, 30 each sliding plane element contains spring and St. V body to model the cohesion and deformation characters, respectively. The St. V body can only slip in the tangential direction. Thus, in the present simulation, a certain percentage of bonded contacts with low cohesion and much bigger tensile strength are installed between two particles to play the role of initial "sliding planes." As the tensile strengths of these "slide planes" are big enough, they could only fail due to the increase of shear stress s and would not fail in tension. In the loading process, it would first fail when the shear stresses exceed the slip precondition (failure criterion). In the present work, different from Wang et al's model, the slip behavior of slide planes is not only related to the cohesion but also the frictional resistance. Thus, after the failure of that slide plane, it would slip when the shear stress exceeds the frictional strength (details are shown in Section 2.1.2) and then the inelastic deformation (ie, slide displacement) appears. The slide displacement would gradually increase with the loading process. When the stress state switched from loading to unloading, as there is a stress interval during which the force of friction changes direction, no displacement occurs then, and the slide plane would stop slipping. 29 Thus, we make a simplification that the slipped slide plane would be rebounded when the unloading starts. The unrecoverable elastic energy would store in the slide planes. As part of the elastic energy stores in the slide plane system, the slip precondition would increase, and the new slip precondition is related to the maximum preloading stress; in other words, the "memory" is generated. In the following loading stage, the plane would slip only when the shear stress exceeds the new slip precondition and cause the change of slope of DRA curve.
Hereafter, there are three basic requirements for the used contact element: (a) It is one kind of bond contact models; (b) it would slip after the slip precondition achieved; (c) it could be rebonded after slipping occurs.
In For the LPBM, the bond body and slip body are in parallel, so the slip behavior is isolated with the bond behavior which could not meet the requirement b. Comparing the LCBM and SJBM with the FJBM, the latter one could set a bigger installation gap, which means it could be locked even big slip displacement occurred after reaching its slide precondition. However, the installation gap for LCBM and SJBM is a fixed value (1e − 6), only when the surface gap is smaller than the mean radius of the particles in its end, it could be bonded. It means that for the LCBM and SJBM, if the slip displacement is relatively big and make the surface gap exceeds the mean radius of the particles, the slipped contact would miss and could not be relocked. In conclusion, only the build-in FJBM could meet the basic requirements to model deformation memory effect in PFC.

| Flat-joint bond model
A brief introduction into FJBM is helpful in understanding the simulation process of deformation memory effect in PFC.

| 3031
TANG eT Al. When the base spherical particle material is generated, the initial flat-joint contacts are installed at the grain-grain contacts according to the relationship between their gaps and the specified installation gap according to Equation (2). When the surface gap g s is smaller than the installation gap g, the contact would be assigned as a FJBM.
where the g is the installation gap, g ratio is the installation gap ratio, and R A and R B are the radius of two contacted particles.
The installed FJBM is composed of several elements, and its characteristics are shown in Figure 4. After a flat-joint contact is installed between two particles, forces and moment are updated based on these elements. The normal force F n e and moment M (e) are computed directly, and the shear force F s e is computed in an incremental fashion. 43 The updating law of the normal force and moment for a flat-joint contact is as follows.
where A is the area of the element, r is moment arm, and is the normal interface stress.
The normal interface stress is given by where k n is the normal stiffness, and g is the interface gap.
When the normal force is updated, it is then used to update (e) in Equation (5).
The normal stress (e) and shear stress (e) acting on an element are given by The shear force-update law for FJBM is as follows: where K s is the tangential stiffness, and s is the tangential displacement. Then, the shear force is used to update (e) in Equation (5).
After updating forces, it is followed by the judgment of crack generation.
The tensile strength is specified as b . If the normal stress b , then the element breaks in tension. As the tensile strength of the initial "slide planes" is set big enough, this law is unnecessary to consider for them.
The procedure considers the shear strength based on the bond state of FJBM as follows.

Unbonded
The frictional shear strength obeys the Coulomb sliding criterion, | ≤ r , the shear stress remains (e) ; otherwise, the shear-strength limit is exceeded, and slip occurs.

Bonded
The shear strength follows the Coulomb criterion with a tension cutoff, where c b is the bond cohesion; and b is the friction angle. If | | (e)| | ≤ r , the shear stress remains (e) ; otherwise, the shearstrength limit will be exceeded, and the bond breaks in shear.
According to the above introduction of FJBM, it could be found that the FJBM is more realistic to model the slide plane than the model proposed by Wang et al because there is an oversimplification for it. 30 The compressive stress normal to the sliding planes has no effect on the slip behavior. The methodology to realize the deformation memory effect in PFC based on FJBM is shown in Figure 5. There is an initial "slide plane" between Ball1 and Ball2 ( Figure 5B). When the shear stress the slide plane suffers exceeds its failure criterion ( (e) > c in − (e) tan b is its initial cohesion value) in the preloading stage ( Figure 5C), the plane reaches its slip precondition. The slip would further occur after the shear stress exceeds its unbonded shear strength (frictional strength) in Equation (7) and results in slip displacement. The slide plane would slip until the end of the preloading stage, and then, it would stop due to self-locking mechanisms (rebonding the slide plane), and residual elastic energy would store in it. The stored elastic energy would increase the c in to a new cohesion value c in + c � ( Figure 5D). Thus, in the following loading process, after the shear stress exceeds its new failure criterion ( (e) > (c in + c � ) − (e) tan b , c ′ is the extra cohesion value results from residual elastic energy) ( Figure 5E), the slip precondition reaches, after which the inelastic deformation would increase again due to the slip behavior. At the start of the first unloading stage, it would relock, which means its slip precondition becomes bigger again ( Figure 5F). As the maximum second loading stress is the same as the first one, it would not slip in the second loading stage ( Figure 5G). Thus, comparing with the first loading stage, there would no inelastic deformation occurring in the second stage. In other words, the slope of the difference strain curve between the first loading and second loading stages would change at the maximum preloading stress point.

| Calibration of microcharacteristics
In the following modeling, LdB granite was chosen as test samples. LdB granite is a typical type of brittle rock, belonging

F I G U R E 5
Methodology and process to model deformation memory effect based on FJBM using PFC3D: (A) Rock model and loading wall; (B) Two ball with one initial "slide plane"; (C) The moment reaching precondition and slip occurs in the preloading stage; (D) At the start of preunloading stage, the slide is relocked and part of elastic energy is stored in it; (E) After reaching precondition, the inelastic strain caused by slide generates in the first loading stage; (F) At the start of first unloading stage, the slide is relocked and its slip precondition becomes bigger; (G) As the second loading stress is same as the first one, the slide plane would not failure to class I brittle rock. Many researchers performed physical experiments on this type of rock. In the present work, the main macro-properties we concerned are Young's modulus, Poisson's ratio, uniaxial compressive strength (UCS), tensile strength, and crack initiation stress. The Young's modulus was obtained by the tangent at the 50% UCS. The crack initiation stress, which indicates the onset of crack growth, is determined to be approximately 40% of UCS by the volumetric strain method for the LdB granite. 44,45 Before DRA method simulation, microparameters calibration is necessary to get a realistic rock model. The uniaxial compression and direct tension test simulations were performed to calibrate the microproperties. In both simulations, the loading velocities of the loading walls are 0.01 m/s, under which the loading condition could be considered to be in a quasi-static state. 46,47 The detail of the calibration procedure can be seen in the cited reference. 48 The microparameters after calibration are listed in Table 1. It should be noted that the crack initiation stress in the PFC simulation is defined as one axial stress value at which a specified number of microcracks induced. This number was defined as 1% of the total number of cracks compared with the failed sample in the present work. 26 After calibration, the simulated macro-properties show good agreement with the laboratory tests (Table 2). Hence, the following DRA method simulations are carried on this calibrated model.

| Simulation for DRA tests
The DRA test simulations are presented in Table 3, and the ratio between the following loading and preloading stress is set as 1.5. To verify the feasibility of this model in the lowstress region, the preloading stress should be lower than the crack initiation stress. The loading was controlled by platen velocity. The loading rate is the same as that in the UCS test simulations (ie, 0.01 m/s). In the simulations, the stop of the unloading stage was controlled by the axial loading value when the axial stress drops to 0.2 times the maximum loading stress.
Besides the high replication of rock materials, another advantage for PFC simulation is the convenience of data measurement. In this simulation, four measurement spheres in the same position of the strain gauges in the experiment were used to measure the strain data. The final strain data were measured using the average strain of these four measurement spheres. 49

| DRA curves with different peak loading stresses
As mentioned previously, based on the established model, the inflection in the DRA curves should be apparent when the test is conducted in the low-stress region where the maximum loading stress is smaller than the crack initiation stress. This behavior is demonstrated in Figure 6, which depicts the strain difference in different directions (axial and lateral directions) vs the axial stress of the established rock model subject to different maximum preloading stress. When the maximum preloading stress is 0.6 MPa, there are no inflection occurs in the DRA curves. This is because the preloading stress is too small, and no slide element reaches its slip precondition. However, there are obvious inflections both in the axial and lateral directions when the preloading stress reaches 2 MPa. In other words, although there would be deformation memory effect when the loading stress is smaller than the crack initiation stress, there still exists a stress threshold below which there would not "memory" generate. Furthermore, when the preloading stress reaches 2 MPa, for the axial direction, the strain differences were first kept at a stable value and then dropped sharply, in which a distinctive "memory kick" can be observed and corresponds to the maximum stress applied in the preloading stage. Along the lateral direction, the differential strains were also first kept at a stable value. The difference observed between axial and lateral directions is that the DRA curve in the lateral direction sharply increases when the applied load reaches the maximum preloading stress. Thus, in the modeling, obvious inflections can be observed. This phenomenon is the same as the experimental results conducted by Hsieh et al. 49,50 Some features of DRA curves would appear when the applied preloading stress is equal or above 3 MPa. Comparing these results with 2 MPa, the differences are as follows: 1. During the loading stage, the strain difference in the axial direction keeps increasing when the stress is smaller than the maximum preloading stress and then sharply drops after it exceeds this value. The strain difference in the lateral direction first decreases with the axial stress and suddenly increases when the loading stress beyond the maximum preloading stress; 2. With the increase of preloading stress, the DRA curves after the inflection become volatile, and the inflection becomes unclear relatively.
The mean values and felicity ratios between the inflection value and maximum preloading stress in different directions are listed in Table 4. In general, although there are some differences between DRA curves, the inflections which represent the "stored-memory" are clear and accurate. It has verified the feasibility of PFC in DRA method modeling when applied with low loading stress.

| Strain energy
As mentioned in Section 2.1, one key point for reproducing deformation memory effect in low-stress region is the lock mechanism of slide planes. The slide element could store part of elastic energy and further cause the increase of slip precondition, which is related to the maximum preloading stress. In PFC3D, the strain energy of the FJBM is obtained by summing the strain energy in each element (Equation 9) which is updated after the force-displacement law 43 : with I (e) = 1 4 r 4 e and r e = √ A (e) where I (e) is the inertia moment of the contact element, and r e is the element effective radius. The k n , k s , I (e) , and A (e) are constant value. Thus, the measured strain energy could reflect the force development of the rock sample system.
The development of strain energy when the maximum preloading stress is 2 MPa is shown in Figure 7. At the end of the unloading stage, there is a strain energy gap between the loading and unloading stage in the pre and first cycles, which means some strain energy was rest and stored in the rock model system (ie, initial "slide planes"). In other words, according to Equation (9), there would be a prestress produced inside the initial planes in the precycle, which is equal to the increase of their slip precondition. The new slip precondition is related to the maximum preloading stress. When the loading stress in the first cycle exceeds the maximum preloading force, the initial "slip planes" would restart to slip and cause some strain energy been stored again. The slip precondition would accordingly increase again. However, as the maximum second loading stress is equal to the first one, the "slide planes" would not slip in the second cycle, and there would be no rest strain energy been stored. Thus, there is no strain energy gap in the second cycle.

| Slip energy
The calculation of the slip energy for the FJBM in PFC is shown in Equation (10). 43 As the (e) c and A (e) in the equation are the constant values, the change of slip energy could represent the development of the slip displacement, that is, the inelastic strain. The deformation memory effect under small This history of the slip energy has been recorded during the whole simulations, and the results are shown in Figure 8. When the maximum preloading stress is 2 MPa, the slip energy first increases when the loading stress increases to approximately 0.75 MPa, and it means that some slide planes reach their slip precondition, and then they slip. At the start of the preunloading stage, those slipped "slide planes" were locked and stop slip. Thus, the slip energy keeps almost unchanged then. In the subsequent loading stage, the feature of the slip energy curve is similar to that of the former one, as it first keeps unchanged and then sharply increases. The slope of the slip energy curve starts to sharply change when the loading stress reaches 2 MPa, and there is a clear inflection corresponding to the pre maximum loading stress. Unlike the preunloading stage, there is also an obvious increase for slip energy in the first unloading stage, which is due to the fact that the planes would slip in the reserve direction in the unloading stage, and they may also reach their slip precondition in this stage. When some slide planes reach their slip precondition in the first unloading stage and further slip, the slip energy will increase.
The overall development law of slip energy for the maximum preloading stress with 3 MPa is similar to 2 MPa, except that the slip energy increases continually in the pre unloading stage. The inflection is still clear in this case (3 MPa). However, with the further increase of maximum preloading stress, the inflection point in the first loading stage almost disappears, which is due to the significant slip behavior of the slide plane in the preunloading stage. When the loading state changes into the unloading state, the previously slipped slide planes would stop slipping due to the self-locking mechanisms, and then the shear force would increase in the reverse direction. During the reverse process, the "memory" inside these planes would lose when reaching their slip precondition in the preunloading process. Those slide planes would not be rebonded, which means they would continuously slip and cannot reach the slip precondition again when the stress in the first cycle exceeds the maximum preloading stress. Thus, there would not be a slip energy inflection point occurs.
The above-mentioned conclusions can also be used to explain the different phenomena between DRA curves with different maximum preloading stresses described in Section 3.1. As the maximum preloading stress is smaller than the crack initiation stress, the inelastic strain is attributed to the frictional slip. Thus, the strain difference obtained from Equation (1) is related to the difference of slip displacement between the first and second loading stage. When the maximum preloading stress is small (2 MPa in this case), no slip behavior occurs in the first loading stage before the loading stress reaches the maximum preloading stress. After reaching the maximum preloading stress, slip behavior occurs, and inelastic strain generates. Although there is slip behavior in the second loading stage, the slip displacement difference of them before the maximum loading stress is smaller than that after the maximum preloading stress. Thus, the strain difference almost keeps unchanged at first. In addition, as the strain difference after the maximum preloading stress is relatively large, the strain fluctuation appears to be slight.
With the maximum preloading stress increases, the reversed slip behavior of "slide planes" also increases both in the preunloading and first unloading stage. It causes a decrease in the slip behavior after the stress exceeds the maximum preloading stress. In other words, some "slide planes" could not store "memory." The inelastic strain (ie, slip displacement) difference before and after the maximum preloading is not large enough to make the strain difference in the initial stage looks unchanged. Generally speaking, for one kind of rock (eg, with certain filler property to control the slip precondition), there exists a suitable range using DRA method for low in situ stress measurement. If the stress is lower than the minimum threshold, no inflection point will appear, while the inflection will become relatively unclear with the increase of stress. However, the accuracy of the DRA results would almost not been affected as long as the stress exceeds the minimum threshold.

| PARAMETRIC STUDY
According to Section 2.1.2, the rock properties obviously affect the DRA results. However, for DRA method, the influence of variations in material properties on the measured results is still unclear. Due to the complexity of rock materials, it is necessary to explore it first based on numerical simulation. For the established model, the key point to realize the deformation memory effect in low-stress region is to model the slip behavior of initial "slide planes." The microparameters that have a major effect on the slide plane property are friction coefficient, cohesion, dip angle of "slide planes," and their fraction. In reality, the cohesion is closely related to the slip precondition, as a larger cohesion leads to a higher slip precondition. The friction coefficient determines when the slide would slip after reaching precondition. The slide plane fraction F represents the initial damage level of the rock sample. The dip angle of initial slide planes is related to the rock anisotropy and its microstructure. Therefore, the parametric study focuses on these four parameters. When performing the parametric study, these parameters are changed one at a time while keeping the other parameter values constant, as shown in Table 5.
For each parameter study, we chose five different values. Therefore, there are 20 models used. The loading strategy and loading velocity are the same as those in Section 2.3. The maximum preloading and first loading stress are set as 2 and 3 MPa, respectively.  Figure 9 shows the DRA curves of rock models with different cohesions in simulations. When the cohesion exceeds 2 MPa, the results are similar to that with 1.5 MPa as there is no inflection corresponding to the in situ stress. Thus, those results with cohesions over 2 MPa are not presented here in order to avoid repetitions. When the cohesion is 0.6 MPa, the strain difference first increases and then sharply decreases in the axial direction. For the cohesion with 1 MPa, the DRA curve first keeps stable and then sharply decreases. This phenomenon is consistent with the experimental results obtained by Attar et al, 11 and their results are shown in Figure 10.

| Effect of the cohesion
Although the maximum preloading stresses are the same (2 MPa), the DRA curves are different: It increases first in Figure 10A while keeps stable in Figure 10B. This phenomenon may be attributed to the different cohesions between two rock samples, which would cause the slide planes inside the sample with small cohesion ( Figure 10A) are easier to slip in the preunloading stage, and further leads to "memory loss." Thus, the DRA curve would first increase while not keeps almost unchanged. The strain difference rate of the DRA curve when the loading stress is lower than the preloading stress is shown in Figure 11. There is a dramatic decrease when the cohesion exceeds 0.6 MPa. When the cohesion value is much lower than the in situ stress, the DRA curve becomes inclined. Besides this, the DRA curve in the axial direction is more easily to be affected. The felicity ratio between maximum preloading stress and the measured value is close for cohesions with 0.6 and 1 MPa. This illustrates that if the in situ stress is in a suitable range (bigger than the minimum threshold value), this technology could be used to measure in situ stress with high accuracy.
However, it could be found that there is also a clear inflection when the axial stress is about 2.25 MPa with the cohesion of 1.5 MPa. Apparently, this inflection is inaccurate. Although this in situ stress data could be used in engineering applications because an overestimated value has a conservative estimation effect on engineering design, it is necessary to understand its mechanism. The number of slide planes which slipped in the loading cycles is counted and shown in Figure 12. It can be found that although there are some planes slip in the preloading stage, the number is quite small and only compose about 0.5% of the total slide planes. It means that these slide planes could not store enough information to generate memory in the next loading cycle. Thus, there would be no inflection corresponding to the maximum preloading stress. It can be found in Figure 12 that the inflection value (2.25 MPa) is corresponding to the point when there are some planes start to slip in the first loading stage. As the number of slip planes is corresponding to the stress value when there are enough slide planes start to slip.
According to the above analysis, there is a question that should be paid more attention when using DRA method. When the DRA method is adopted to measure the in situ stress, the actual in situ stress could make the slide planes reach their slip precondition. However, if the in situ stress is not large enough to make adequate planes to slip, it would produce overestimated measurement results. In order to judge the veracity of the measured value in the engineering application for determining the in situ stress in low-stress region, the following checking procedure is needed:  Figure 13 presents the effect of the friction coefficient of "slide planes" on DRA curves. The results with friction coefficients of 0.4, 0.6, and 0.8 are similar to that of 0.2: DRA curves first keep stable and then drop or ascend sharply when the axial loading stress reaches the maximum preloading stress. Thus, in order to avoid repetitions, similar results are not given. The difference between DRA curves with friction coefficients of 0 and 0.2 is that for friction coefficient of 0, the strain difference first increases while for friction coefficient of 0.2, it keeps stable. The strain difference rate is shown in Figure 14. Except for the case with a friction coefficient of 0.0, the rate is almost near 0, which means the strain difference curve is almost invariable before reaching preloading stress. Same as the cohesion parameter, the axial direction result is easier to be affected. The friction coefficient parameter would affect the slip behavior after the planes reach their slip precondition. Thus, the shear stress of a certain slide plane (id = 192) with different friction coefficients was recorded in the precycle and shown in Figure 15. The solid and dotted lines represent the loading and unloading stages, respectively. After the slide plane reaches precondition, the shear stress would sharply decrease to a value which is related to the friction coefficient and large the friction coefficient is, large this value is. Then, the shear stress would increase except the case with a friction coefficient of 0. Thus, the shear stresses increase in the order 0.8, 0.6, 0.4, 0.2, and 0.0. The large the shear stress is in the initial of the unloading stage, the more difficult for the plane to slip in the unloading stage. For the model with a friction coefficient of 0, the shear stress of the slide plane keeps at 0 from it reaching slip precondition until the end of the loading stage. Thus, the slide plane would be easier to slip in the unloading stage and lose the stored memory. In the following loading stage, this slide plane would not reach its slip precondition again. The slip displacement difference before and after the maximum preloading stress is not large enough to make the strain difference in the initial stage looks unchanged. For the models with friction coefficients larger than 0, the strain differences keep unchanged in the initial stage. However, although there are some differences, the accuracy is high. Thus, the friction coefficient has little influence on the accuracy of the DRA method.

| Effect of dip angle
In this section, the effect of the dip angle of slide planes on the DRA results will be analyzed. The dip angle in the PFC simulation represents the angle between the normal direction of the slide plane and the z-axis. The dip angle of all the slide planes was set at a certain range (Table 5) in these five cases, as shown in Figure 16. The law of DRA curves of Case b-d is almost the same: The strain difference first keeps invariable at a fixed value and then changes sharply accompanied by a clear inflection. Thus, only the DRA curves of Case a, b, and e are shown in Figure 17.
The dip angle has a significant effect on the DRA results. When the slide planes are vertical to the loading direction (Case a), there is no inflection; when the slide planes are parallel to the loading direction (Case e), the inflection is not clear; in other cases (Case b-d), the inflection is clear, and the measured value is accurate. The strain difference rate is shown in Figure 18. The rate first increases and then decreases with the dip angle; however, the value for these cases is small. The rates in the lateral direction are closer to 0 compared with the ones in the axial direction.
As rocks can be regarded as one kind of granular mediums 35 composed of particles which displace independently from one another and interact only at contact points, the interaction between two particles in the axial loading condition is simply shown in Figure 19. Thus, the shear force that the slide plane subjected is: where is the dip angle of the slide plane. Thus, As the dip angle is small in Case a, the shear stress at the slide plane is too small to cause it to slip. Thus, there is no inflection point. For Case e, due to the big dip angle, the shear stress is too large and can cause some relocked planes to slip in the unloading stage. These slide planes would lose "memory" and make the inflection unclear. However, besides these two cases, we could get accurate in situ stress data. According to the above analysis, when the tested rocks are layered or stratified rocks, it is suggested not to drill the rock sample from rock cores in the direction parallel or vertical to the layer direction. Figure 20 shows the DRA results of rock samples with different fractions of "slide planes." The fraction represents the initial damage of rock samples. There is a clear inflection in DRA curve when the fraction is greater than 0, and the accuracy is high (similar results with the fraction of 0.1, 0.15, and 0.2 are not given for avoiding repetitions.). The strain difference rate is shown in Figure 21. The distribution of the strain difference rate in different cases is similar except when the fraction is 0%. The law in different directions is the same as the previous ones. Combined with the results in the Section 4.1, if the fraction of slipped slide planes is much low, no inflection will occur. Thus, if there exists rock with a little slide planes (eg, weak plane), it is not capable of storing enough "low-stress memory."

| CONCLUSION
In the present study, a method to model the deformation memory effect of rock in the low-stress region is proposed within a frame of particle-based discrete element method using FJBM. The pre-existing bonded contacts inside the rock sample with large tensile strength and relatively small cohesion are created to model the cohesion property of initial "slide planes." Besides the cohesion property, the friction slip behavior and self-lock mechanism are introduced to realistically model the "slide plane." The simulation results (ie, DRA curves) have proven that the DRA method in low-stress region is valid in the numerical modeling environment. The following conclusions can be drawn: 1. There exists a suitable range of in situ stress in low-stress region for DRA method. Within this range, there is a clear inflection corresponding to the maximum preloading stress, and the accuracy is high both in the axial and lateral directions. Although there would be deformation memory effect when the loading stress is smaller than the crack initiation stress, there still exists a stress threshold value below which the "memory" would not be generated. The threshold is mainly related to the cohesion and dip angle of "slide planes." 3. If the in situ stress is much larger than the stress threshold (still smaller than the crack initiation stress), the inflection will become unclear with the increase of in situ stress due to the reslip behavior in the unloading stage which would make part of "memory" lost. 4. Four factors influencing the deformation memory effect are discussed, including cohesion, friction coefficient, dip angle, and fraction of slide planes. It is found that the cohesion and dip angle are highly related to the threshold value. For intact rock (with little fracture), it is not capable of storing "low-stress memory." If there exists enough fracture inside the rock sample, the fracture level would not influence too much on the DRA test results. In addition, for the friction coefficient parameter, it has a slight influence on the deformation memory effect. 5. Even there may exist a clear inflection, the measured value might not be right if the in situ stress does not exceed the minimum threshold value. Thus, when using DRA method in low-stress region in the engineering application, other measurement methods should be firstly adopted to confirm its feasibility. 6. It is suggested not to drill the rock sample from rock cores in the direction parallel or vertical to the layer direction when the tested rocks are layered or stratified rocks.