Journal of Geophysical Research: Solid Earth Quantiﬁcation of Fracture Interaction Using Stress Intensity Factor Variation Maps

Accurate and ﬂexible models of fracture interaction are sought after in the ﬁelds of mechanics and geology. Stress intensity factors (SIFs) quantify the energy concentrated at the fracture tips and are perturbed from their isolated values when two fractures are close to one another. Using a three-dimensional ﬁnite element fracture mechanics code to simulate static fractures in tension and compression, interaction eﬀects are examined. SIF perturbations are characterized by introducing three interaction measures: the circumferential and maximum SIF perturbation provide the “magnitude” of the eﬀect of interaction, and the ampliﬁcation to shielding ratio quantiﬁes the balance between increased and decreased SIFs along the tip. These measures are used to demonstrate the change in interaction with fracture separation and to ﬁnd the separation at which interaction becomes negligible. Interaction maps are constructed by plotting the values of the interaction measures for a static fracture as a second fracture is moved around it. These maps are presented for several common fracture orientations in tension. They explore interaction by highlighting regions in which growth is more likely to occur and where fractures will grow into nonplanar geometries. Interaction maps can be applied to fracture networks with multiple discontinuities to analyze the eﬀect of geometric variations on fracture interaction.


Introduction
Fracture interaction in brittle rocks is of great interest at a variety of scales. Fractures alter the surrounding stress state of the medium, enhancing or reducing stress depending on the mode of deformation, location, and fracture orientation (Olson & Pollard, 1991). Predicting the behavior of a fracture system requires a thorough understanding of the stress state within the material and the current fracture geometries. This is not straightforward to characterize in many settings, particularly for geological materials due to their heterogeneity, anisotropy, and the uncertainty of subsurface characteristics relative to exposed outcrops. At smaller scales, linkage of fractures has a large control on the hydraulic properties of rocks and reservoirs, with implications for hydrocarbon production (Nelson, 2001) and radioactive waste disposal (Tsang et al., 2015). At larger scales, interaction effects are a key control over the formation of geological structures (Delaney & Pollard, 1981), rift evolution (Gawthorpe & Leeder, 2000), and fault slip during earthquakes (Green et al., 2015;Zachariasen & Sieh, 1995). Determining when faults transition from soft to hard linkage requires accurate observations and a detailed history of displacement (Duffy et al., 2015).
When a fracture is present in a tensile stress field, it alters the stress in its local area in two ways (Kachanov, 1993;Price & Cosgrove, 1990;Rives et al., 1992). First, a fracture generates a shielding zone around its surface, also known as a stress shadow, in which local stress is reduced. Second, at the tip of a fracture, an amplification zone is generated, also known as the concentration zone, in which local stress is increased. When another fracture enters these zones, this change in stress results in significant changes to its behavior. When the fractures are similar in size, they induce a reciprocal effect upon one another. In full 3-D space, this interaction is complex, and the change in stress varies significantly over each fracture, particularly when fractures are overlapped. Interaction is important because it may both enhance and preclude growth and is a controlling mechanism of the self-organization of fracture patterns. In turn, self-organization may influence the ensuing mechanical and flow properties of the fracture pattern (Olson, 2004;Olson & Pollard, 1991).
Numerical simulations are instrumental in understanding fracture interaction. A variety of techniques have been utilized to numerically model fractured geological media (Jing & Hudson, 2002), many of which include discrete fractures. Three main approaches can be distinguished: stress based, debonding based, and stress intensity factor (SIF) based. Stress-based methods usually rely on displacement fields generated by the finite element method (FEM) (Carter et al., 2000;Pouya, 2015;Rabczuk, 2013;Riahi et al., 2010). They rely on measurements of stresses local to the fracture tips and therefore require significantly refined meshes. Debonding-based methods are routinely used in conjunction with the distinct or discrete element methods (Lei et al., 2017;Lisjak & Grasselli, 2014) and model rock mass interaction as agglomerate rock matrices made up of bonded particles, which mimic the development of fractures by releasing element bonds when local stresses are exceeded. Debonding methods follow the preexisting discretization and require the definition of microscale properties using calibration tests that are rock type and scale dependent.
Stress intensity factor (SIF)-based fracture modeling has received great attention in material design to predict failure (Busfield et al., 2005;Kamaya, 2008a;Moussa, 2002;Stonesifer et al., 1993). SIFs provide an effective method for analyzing fracture propagation and interaction (Lam & Phua, 1991;Paris & Erdogan, 1963). They are one of several methods of fracture analysis, quantifying the energy released at a fracture tip during crack growth (Anderson & Anderson, 2005;Irwin, 1956). They can be used to evaluate failure criteria for modeling fracture propagation (Brace, 1960;Lawn, 1993;Paris & Erdogan, 1963). They have become a key part of fitness-for-purpose analyses when considering the impact of interacting cracks on a material's lifetime (Wiesner et al., 2000). SIFs have also been used in geological models to predict fracture behavior (Olson, 2004;Paluszny & Zimmerman, 2013;Philip et al., 2005;Renshaw & Pollard, 1994) but have not been routinely adopted for investigating geological fracture interaction. Kachanov (1987) presented an approximate method for calculating how the SIFs of one fracture are affected by another fracture. The system is simplified by using a transmission factor which depends on the fracture geometry and orientation of the stress field but not the magnitude of the stresses. The 3-D formulation of this method incorporates the idea that fracture interaction may amplify (increase) or shield (decrease) stress locally depending on the geometry. Other early analytical approaches include that of Fabrikant (1987), who used an iterative procedure to evaluate each fracture's effect on one another. More recently, work has expanded on these methods to improve accuracy or investigate more complex systems of fractures. Laures and Kachanov (1991) extended this analysis to 3-D microcrack arrays using the same method as Kachanov (1987), focusing on cases where one fracture interacts with a group of smaller fractures. They show that 3-D interaction effects are generally weaker than 2-D interaction effects. Gross (1982) first introduced the method of polynomial expansions of tractions on fractures, providing approximate SIFs for multiple straight cracks. Benveniste et al. (1989) used a superposition scheme for which Kachanov (1987)'s method is a special case. Zhan and Wang (2006) used a first-order average traction method to study the same set of geometries as Kachanov (1987).
Numerical methods for calculating SIFs provide a convenient means to study SIF alteration, because any arbitrary fracture geometry can be investigated. In 2-D, Fu et al. (2012) demonstrated amplification and shielding using the displacement correlation method to find SIFs. Yan (2010) investigated the interaction of circular arc cracks in tension. Fan et al. (2016) investigated the case of two fractures in compression using the FEM and compared their numerical results to photoelastic and uniaxial compressive experiments. Legrand and Lazarus (2015) used a finite perturbation method to study fracture interaction in 2-D. Along with SIF perturbations, they studied fracture growth and coalescence, evaluating growth using a fatigue law to govern fracture propagation. Kamaya (2008b) used the FEM to study SIF variation along the tips of surface cracks, using the SIF perturbation to accurately predict growth (Kamaya, 2008a). Most of these studies investigate interaction of static geometries and yield interaction relationships that depend either on distance or on relative spatial and geometric arrangements but are restricted to specific geometries such as coplanar or stacked configurations. To aid in interpretation of fracture interaction and growth, there is need for methods that can simulate interaction in full 3-D space, without constraining the size, orientation, or number of fractures. A flexible numerical approach can be readily applied to any fracture geometry, subjected to arbitrary boundary conditions, while analytical and semianalytical approaches may be constrained to specific geometries and conditions.
With the development of advanced 3-D FEM based models for studying fracture propagation, SIF-based analysis can be applied to fracture interaction in order to understand the importance of relative fracture spacing. In the present work, a fracture mechanics-based 3-D FEM model, outlined by Zimmerman (2011, 2013), is used to investigate fracture interaction. Advances in methods for FEM SIF computation provide highly accurate results for closely interacting fractures, in particular, by solving for the interaction integral using a disk-shaped domain integral method (Nejati et al., 2015a), without requiring any special mesh structure around the fracture tips. The present work extends the analysis of fracture interaction by investigating changes in the SIFs for all modes of fracture tip deformation, namely, modes I, II, and III. Characterization of interaction is captured by three novel measures for quantifying the spatial distribution of fracture growth energy. Interaction is quantified based on the perturbation of the SIFs by comparing interacting fractures to undisturbed single fractures in the same stress field. The spatial variation of these interaction measures are proposed as "interaction maps" or "stress intensity factor variation maps," used interchangeably in this manuscript. These maps provide a new way to assess the extent and magnitude of fracture interaction.

Numerical Method
The fracture modeling workflow is summarized in four steps: geometry, meshing, deformation, and SIF computation. In this work, propagation of fractures is not explored; instead, the focus is on understanding how interaction between fractures enhances or suppresses growth as a function of spatial arrangement. The finite element-based fracture mechanics framework used herein is previously detailed by Paluszny and Zimmerman (2013) and has been validated for both 2-D and 3-D geometries (Paluszny & Matthäi, 2009;Paluszny & Zimmerman, 2011). The fracture geometry is stored independently of the mesh, and so the mesh is merely an instrument that is used and disposed at each step, created to satisfy both geometric and solution field constraints of the particular simulation stage.
Fractures are represented as smooth surfaces growing in response to tension or compression, governed by the mechanics of brittle failure. The volume is discretized using a volumetric mesh. An initial fracture geometry is created inside a volumetric box representing the rock mass. Fractures can be assigned any surface shape, including nonplanar geometries. In this work, fractures are assumed to be disk shaped to simplify analysis, but all methods proposed are also effective for fractures of any shape, provided the shape can be meshed with a consistent refinement along its tip. At each step, an "ideal" mesh is constructed for the given geometry. The mesh is composed of two element types: isoparametric quadratic tetrahedra in the volume and isoparametric quadratic triangles on the fracture surface elements. Triangle and tetrahedral elements around the crack tip are taken to be quarter point elements (Banks-Sills & Sherman, 1986), which capture the singularity ensuing at the fracture tip by shifting the midside node one quarter toward the fracture tip. These have been found to significantly improve the measured stress intensity factor values in unstructured tetrahedral meshes (Nejati et al., 2015b). In addition, the created mesh is refined near the fracture tip to improve sampling of the displacement field. This refinement is a function of the geometric discretization of the fracture tip. Automatic remeshing enables many different fracture geometries to be tested without human interaction.
The FEM approximates the deformation field numerically, as a function of stresses, material properties, and geometry, and ensuing stress and strains are derived assuming that the matrix is homogeneous, isotropic, and linear elastic. Stress and strain are governed by where is the Cauchy stress tensor, is the infinitesimal strain tensor, 0 and 0 are the initial stress and strain, and D is a linear elastic stiffness matrix (Cook et al., 2007). For a static system,  Rice, 1968). Using the I integral is beneficial as it can be decomposed into individual SIFs for mixed mode problems (Nakamura & Parks, 1988;Yau et al., 1980). The I integral is widely considered one of the most accurate methods for calculating SIFs in 3-D (Banks-Sills, 2010; Bremberg & Faleskog, 2015;Walters et al., 2005). The method presents a number of challenges that make it difficult to perform on a standard FEM mesh near the crack tip, such as the challenge of integrating over an unstructured mesh and conversion to a curvilinear coordinate system. These challenges can be mostly overcome by resampling the crack tip region and using a disk-shaped domain integral formulation of the I integral, outlined in full by Nejati et al. (2015a). Integral methods have been shown to be valid when the domain includes material changes or discontinuities (Yu et al., 2009), so the disk-shaped domain integral is suitable for study of very close fractures.
This methodology is implemented as part of the Imperial College Geomechanics Toolkit (Paluszny & Zimmerman, 2011), which is built upon the Complex Systems Modelling Platform (CSMP++) API (Geiger et al., 2001) and is mainly oriented at computational rock mechanics and the simulation of fracture growth (Paluszny & Zimmerman, 2013;Salimzadeh et al., 2017). The FEM inversion to calculate equations (1) and (2) is performed by the Fraunhofer SAMG solver (Stüben, 2001), and meshes are generated using a commercial octree-based 3-D package.

Validation of Domain Integral Method
The analytical solutions for the SIFs along the tip of an inclined penny-shaped fracture under uniform tension are provided by Kassir and Sih (1975): where K is the stress intensity factor (SIF) in units of Pa m −1/2 ; I, II, and III are the crack tip deformation modes; is the applied normal far-field stress; r is the fracture radius; is the angle around the fracture; is the inclination of the fracture plane relative to the load; and is the Poisson's ratio of the medium. K I is constant around penny-shaped fractures of any inclination and equal to 0 when the fracture is parallel to the direction of stress. K II and K III change sign on either side of the fracture (top and bottom for K II and sides for K III ). Stress intensity factor error, % Disk with quarter points Disk without quarter points DC with quarter points DC without quarter points Figure 2. Error in stress intensity factors as mesh size changes. The error, calculated using equation (6), is shown for a fracture with the same geometry as Figure 1 for different numbers of fracture tip nodes. The number of elements in the mesh (from 8,561 to 69,816) scales roughly linearly with the number of tip nodes (from 8 to 160) when there is one fracture in the volume. E = 1 Pa m −1/2 , and normal stress = 1 Pa. The results match the analytical solution very well. Small scatter is found on all the modes as a result of the mesh discretization.
where i is one of the three deformation modes (I, II, or III), L f is the fracture circumference, K A i is the analytical solution for the mode i SIF, K N i is the numerical solution for the mode i SIF, and dl is the distance between fracture tip nodes. Figure 2 shows the variation in e against the number of tip nodes on the fracture, which is used to control mesh refinement. The number of elements in the volume scales linearly with the number of tip nodes. The disk-shaped domain integral method ("disk") is compared to the displacement correlation ("DC") method, another SIF computation method which is more common in the fracture mechanics literature (Branco et al., 2015;Shih et al., 1976). The results with and without quarter points demonstrate how quarter point elements improve accuracy. The integration radius is set to half the tip length on the fracture, as this produces an accurate result at a good range of refinements. Fracture tip nodes are evenly spaced, so the domain radius can be calculated from the fracture radius.
For higher density meshes, the disk-shaped domain integral technique is shown to be more accurate: the smallest error of 0.9% is found at 120 tips. Eighty tips was considered the best balance between accuracy and computation time, with an error of 1%.

Stress Intensity Factor Changes Due To Interaction
The manner in which a pair of fractures interacts is investigated by characterizing how SIFs change when two fractures are close. These changes depend on fracture proximity, size, relative position, and the type of stress applied and have been widely observed (Fan et al., 2016;Kachanov, 1987;Kamaya, 2008b;Legrand & Lazarus, 2015;Olson & Pollard, 1991;Yan, 2010). The key method of examining interaction, explored by other authors (Kamaya, 2008b;Legrand & Lazarus, 2015;Yu et al., 2009), is to first calculate the SIFs around a fracture that is interacting with others (K int i ). Then, all but one of the fractures from the domain are removed, and SIFs are calculated for that isolated fracture (K iso i ). In this way, the same geometry for one fracture is tested twice, and the changes in the SIFs can be examined through the ratio between K int i and K iso i . The two different fracture domains will have different meshes, but their tip nodes are in the same locations for both calculations, and the mesh around each fracture has a similar structure. The fracture spacing convention introduced by Kachanov and Laures (1989) is used in this work, where the spacing is scaled by the fracture radius and reported in the form Δ∕2r.
In this work, it is assumed that any difference between K int i and K iso i is due to interaction effects. However, both K int i and K iso i will be subject to discretization errors, creating additional small differences. By using a numerical result with the same local mesh around the fracture, similar magnitude discretization errors are present at each tip in both K int i and K iso i . Therefore, the ratio K int i and K iso i calculated at each pair of tips factors out these errors. While replacing K iso i with an analytical solution reduces the number of terms which contain numerical error, this would compound discretization errors along with SIF perturbations. Additionally, analytical solutions for isolated SIFs are only available for specific geometries and boundary conditions, whereas fractures can be simulated numerically for any orientation.

Interaction Compared to Other Methods
No exact analytical solution is available for the change in SIFs due to interaction. Instead, results can be compared to the perturbations found through other semianalytical or numerical methods (Legrand & Lazarus, 2015). The first method for comparison is Kachanov and Laures (1989), whose results provide an approximate  Kachanov and Laures (1989) and Legrand and Lazarus (2015) are provided for the same geometry. (b) The same data from Figure 3a, but with the y axis limits changed to (1.0-1.1). Lines have been placed to connect Kachanov and Laures's (1989) points for ease of visualization.
solution that closely matches other semianalytical methods for calculating SIFs for interacting fractures (Fabrikant, 1987;Gross, 1982). The second method utilizes the results of Legrand and Lazarus (2015), obtained using the finite perturbation method. Figure 3 shows the SIF perturbation quantified through the ratio K int i ∕K A i , obtained using these two methods, and the present method, for a fracture under tension when a second coplanar fracture is nearby. The fractures are oriented perpendicular to the stress field and are shown for two separations, Δ∕2r = 0.05 and Δ∕2r = 0.00025. Both fractures have 160 tip nodes. At = 0 ∘ , where the two fractures are closest together, K int i is higher than K A i as a result of interaction. At the Δ∕2r = 0.05 separation, K int i ∕K A i = 1.3 and decreases as the angle increases around the fracture. At separation Δ∕2r = 0.00025, K int i is significantly higher, and K int i ∕K A i = 3.2. K int i decreases at a much faster rate with angle at the close separation. At = 60 ∘ , both separations converge to similar values of K int i ∕K A i . At 160 ∘ , the far side of the fracture is almost unchanged by interaction, yielding values very close to the isolated solution.
All three methods yield a very similar trend. The present simulation results match those of Kachanov and Laures (1989) particularly well, except for the low separation case where the fractures are very close ( < 20 ∘ ). In this region, the present simulation results have higher K int i ∕K A i and match Legrand and Lazarus's (2015) results closely. Δ∕2r = 0.00025 is a particularly challenging case to discretize, because the region between the . Example mesh around two fractures with separation Δ∕2r = 0.05. Each fracture has 30 tip nodes. Volumetric elements are cut by the visualization plane, so some appear to be thin and elongated, whereas they actually pertain to a well-formed tetrahedral element.
fractures becomes very small relative to the fracture size. To place sufficient elements in this region requires an unreasonably high mesh density. Discretized methods that rely on sampling of this region produce scattered SIF values where the fractures are close, which can be seen in Figure 3 for the Δ∕2r = 0.00025 case. Similar scatter was reported by Legrand and Lazarus (2015) for Δ∕2r < 10 −4 . Figure 4 shows an example mesh used in Figure 3 with 30 tips and Δ∕2r = 0.05, with two fractures in the model to demonstrate the way refinement is focused on the fracture tip. Figure 5 shows K int I ∕K iso I for six fracture pairs in two geometries to demonstrate how the magnitude and sign of the interaction changes as the fracture orientation changes. Figures 5a-5c show coplanar interaction where mode I is amplified along the closest part of the fracture. As the separation between fractures increases, K int i ∕K iso i is much lower, becoming effectively equal to unity at Δ∕2r = 1.50. Figures 5d-5f show a stacked configuration where mode I is uniformly shielded as a result of interaction. This effect is diminished with increasing separation but could be said to have a larger range than the coplanar geometry, because, at Δ∕2r = 1.50, mode I is still reduced, as observed by other analyses of stacked cases (Kachanov, 1993). In the stacked configuration, mode II varies in a similar but opposite manner to mode I, where it is amplified rather than shielded. Mode II is perturbed uniformly all the way around the fracture, because the fracture tips are all the same distance from the interacting fracture.

SIF Perturbations for Common Geometries
Coplanar and stacked geometries are the "end-member" cases for fracture interaction, where pure amplification and pure shielding take place, respectively. The transition between these cases is the overlapped case, where part of the fracture is shielded and part is amplified. As the fractures move from coplanar to stacked, the mode I SIF change becomes negative at the point where the fractures are overlapped and positive where they do not overlap. Mode II and III SIFs also change their perturbation behavior across this transition, giving rise to the changes in growth angle which are well known to occur as fracture offset changes (Thomas & Pollard, 1993). This behavior, while intuitive for many applications of fracture mechanics, can be quantified and better understood through further analysis of SIF perturbations. Moreover, the relative orientation of the fractures will further affect the manner in which interaction takes place (see section 5).

Measures of Interaction
In this section, several measures are proposed to quantify and characterize interaction. These "measures of interaction" derive from simple relationships similar to K int i ∕K iso i . While the ratio K int i ∕K iso i is powerful on its own, it cannot be expanded to modes II and III for tension, as these modes change sign around the fracture, for example, Figure 1. Therefore, use of K iso i in the denominator should be avoided, as it approaches 0 for certain angles. The same is true for mode I in certain fracture orientations or stress states. Additionally, K int i ∕K iso i is also only measured on individual tips, with many values per fracture.
It would be useful to have a single quantity that captures how a fracture is changed due to interaction. Summarizing interaction with one value enables analysis using interaction variation maps, aiding understanding of past and potentially arising fracture geometries. For simulations containing hundreds and thousands of fractures, cheap and local measures of interaction will aid the discretization and level of detail required at different locations of the mesh. The methods proposed herein for characterizing interaction are effective for all deformation modes and boundary conditions. Three such methods are proposed as follows. The first, "circumferential SIF perturbation," computes the total area between the K int i and K iso i curves, either as as the sum of all modes (C) or for one specific mode (C i ). The second, "maximum SIF perturbation" (M i ), computes the maximum perturbation of K i for the fracture. Both of these methods are used to quantify the magnitude of the interaction. The third, , the "amplification to This measure can be related to one specific location on the fracture tip at which the difference was largest. This makes M i a local measure of fracture interaction. The maximum SIF perturbation is unique between the three deformation modes, as the tip node with the maximum M I is not necessarily the same as the tip node with maximum M II or M III . The magnitude of M i is tied to . As the method relies on the value of one particular tip, rather than averaging over a group of tips, it may be subject to greater discretization errors.
An illustration of the differences between C and M i is shown in Figure 6a. This plots the mode I SIF as a function of the tip angle in the same manner as Figure 1. The shaded region shows the area between the interacting and isolated curves (C I ), and M I is shown as the maximum change. Both C and M i use the absolute values of K int i and K iso i in their calculations. Thus, they only provide information on the magnitude of interaction and not whether it is amplifying or shielding.

Definition of the Amplification to Shielding Ratio, i
When measuring interaction, it is important to consider whether stress has been amplified or shielded. In the simple cases of coplanar interaction, stress is only amplified, for example, Figures 5a-5c. In cases where the fractures are stacked and at least partially overlap, stress is shielded, for example, Figures 5d-5f. The type of interaction, either shielding, amplification, or a mixture of the two, is quantified by the amplification to shielding ratio ( i ), which is formulated as follows. Two lists of values, a i and s i , are created from the difference between the interacting and isolated SIFs. a i contains the values where the SIFs for mode i have been increased, and 0 otherwise, and s i contains values where the SIFs for mode i have been decreased, and 0 otherwise, THOMAS ET AL. These lists therefore contain information for the tips which have been only amplified ( a i ) and only shielded ( s i ). The areas under the curves of a i and s i are then calculated:

QUANTIFICATION OF FRACTURE INTERACTION
The balance between amplification and shielding can be found by comparing these two areas. Thus, i is defined as Therefore, i equals 1 when mode i is purely amplified and equals −1 when purely shielded. When i > 0, the SIFs at most of the tips of the fracture are amplified by the presence of another fracture, and when i < 0, most of the tips are shielded.

Comments on Interaction Measures
Numerical simulations inevitably produce different results when changes are made in the domain, even if those changes are slight, due to discretization and numerical errors, as discussed in section 3. This means that there will always be differences between K iso i and K int i , even when the meshes are identical and fractures are very far apart. Therefore, the measures of magnitude (C and M) do not converge to 0 when fractures do not interact but do become very small. This is not a problem for analysis of C and M i , because the size of this error is small compared to cases when interaction is significant. However, scatter between the values of K iso i and K int i will result in i being essentially random. When analyzing a set of i values in space, a threshold should be used to ignore the value of i if either C i or M i suggests that there is very little interaction. Conveniently, C i has the same magnitude regardless of the SIF magnitude because it compares the ratio of two areas. A threshold of C i = 0.005 is used in this work.
The proposed interactivity measures share some commonalities with the transmission factor used as part of Kachanov's (1987) analytical method. For coplanar fractures, the transmission factor is increased by the presence of the other fracture. In 3-D, the sign of the transmission factor changes to correspond to cases of fracture interaction either amplifying or shielding stress depending on the geometry. For fracture arrays, the transmission factor concept is extended to an interaction matrix, containing transmission factors which describe the effect of each fracture on one another. The interaction measures suggested here are similar but are now separated into magnitude and sign of interaction, where C or M i can be used to show the magnitude of interaction, and i describes the balance of amplification or shielding. These measures are also assigned one per fracture, rather than one per fracture pair.  Figure 7 shows how C and M I vary in the simple case of increasing fracture separation for three geometries. C is used rather than C I to highlight the effect of combining the three modes together. All graphs (Figures 7a-7c) use the same axes, with C on the left and M I on the right. In all three geometries, C and M I , demonstrate how interaction reduces with separation. Both interaction measures show that shielding interaction when fractures are stacked ( Figure 7c) has a larger range than coplanar geometries (Figures 7a and 7b). The separation at which C falls below 0.005 is marked with a red dashed line, showing that shielding effects have a much longer range than amplification effects.

Variation of C and M With Separation
For the same separation variation of two fractures, their interaction is shown to not only depend on their spacing but also on their relative positioning with respect to the far-field stresses and their relative positioning with respect to one another. The difference in the behavior of C and M I in each case is of particular interest, especially when comparing the two coplanar cases. In Figure 7a, only mode I is changed by interaction, because the fractures are perpendicular to the stress field and K II and K III are 0. Figure 7b has the same geometry, but is oriented at 45 ∘ to the stress, meaning modes II and III are nonzero and are perturbed by the interaction. M I only relies on mode I, so is higher in Figure 7a where the interaction is confined to one mode and smaller in Figure 7b), when the interaction is shared by three modes. C is the sum of the SIF perturbation for all three modes, so it has a higher magnitude in Figure 7b than in Figure 7a. Discretization errors cause small amounts of scatter and are much higher in M I than C.
In summary, C and M i are both capable of capturing the decrease in interaction with separation; however, care should be taken to examine which modes are being altered in each case. Separating C into its mode-specific values C I , C II , and C III to summarize the change in each mode, using equation (8), may be preferable during interaction analysis, in order to identify where fractures will become nonplanar during growth (Thomas & Pollard, 1993).
For numerically derived values to be useful for analysis, they must not vary strongly as a function of mesh refinement. Figure 8 shows

Interaction Maps
In this section, fracture interaction maps are constructed to analyze the effect of one fracture on another, specifically by capturing SIF variations comprising multiple interaction cases, within a single image. They show the variation of C i , M i , and i as the relative position between two fractures changes. A fracture interaction map is constructed by performing simulations in which a base fracture is placed in a fixed position, perpendicular to the direction of stress. A second fracture of equal radius is placed nearby with the same orientation. The values of C I-III , M I-III , and I-III are then measured on the base fracture in tension. The second fracture is then moved to other locations, and the same measurements are made. This creates a grid of values in the x, z plane around the fracture, which is interpolated to create a continuous map of SIF variation around a static fracture. Figure 9 demonstrates this process for three locations alongside the fracture interaction map for I .
The objective of creating an interaction variation map is to explore the variation of the magnitude and type of interaction, that is, amplification, shielding, or a mix of the two, depending on relative fracture position and orientation. Compared to plotting the change in stress around a fracture (e.g., Grechka and Kachanov, 2006), interaction maps observe the change in the SIFs. Therefore, the impact of interaction is more directly quantified in terms of changes at the fracture tip and can be examined separately in modes I, II, and III.  This process could be used to create interaction maps of any interacting fracture orientation or size, as long as the sampled locations do not cause fractures to intersect.
Fracture interaction variation maps are presented in Figures 10-12 showing values of C I-III , M I-III , and I-III , respectively. Each figure includes maps for two orientations of the interacting fracture. Figures 10a, 11a, and 12a show the interaction map generated from two fractures inclined perpendicular to the direction of tension, and Figures 10b, 11b, and 12b show the map with the interacting fracture inclined 45 ∘ toward the other fracture. Figure 10 shows fracture interaction maps of C I , C II , and C III for two geometries. The coplanar geometry (Figure 10a) shows the gradual change between coplanar and stacked geometries. The magnitude of the interaction in the stacked configuration, along x = 0 as shown in Figure 7c, was shown to be higher than in the coplanar configuration, along z = 0, as shown in Figures 7a and 7c. This behavior is clear on the interaction variation map of C I , where the extent of the highly perturbed region is much larger when the fractures overlap and enter one another's stress shadow.

Interaction Variation Maps of C i
In Figure 10a, C I shows that mode I is the most perturbed interaction mode, particularly when fractures are overlapping. Modes II and III are perturbed at approximately half the magnitude of mode I. Mode II perturbations are not increased uniformly in the overlap region. This is because mode II tends to be perturbed only where fracture tips are close to one another. Otherwise, mode II relatively unperturbed if the tip lies above the fracture surface. This leads to larger magnitudes when more of the fracture tips are lined up vertically, creating a region of lower C II between x = 0.5 and 1. These changes in C II highlight where fractures will undergo nonplanar growth as a result of interaction, due to increased mangitudes K II (Cotterell & Rice, 1980), as observed in en echelon fracture sets (Delaney & Pollard, 1981). Mode III is most perturbed when the fractures are partially overlapped and is unperturbed when the tips are close together. This results in its highest value being in the region of x = 1, z = 0.5. Outside of this overlapped region, C III is essentially unperturbed; that is, it falls below the 0.005 threshold discussed in section 4.1.4.
The inclined geometry in Figure 10b shows much lower magnitudes of C. This is partially because the fracture surfaces cannot be placed exactly next to one another, so the highly interacting stacked case cannot be produced at any interacting fracture location. Regions where mode I is more strongly perturbed are present when the interacting fracture's tip is very close to the base fracture plane. In the same manner as in the coplanar case, mode II is larger when the fracture tips are close and overlap in Figure 10b, creating a trough along the line z = 0.75 for the 45 ∘ inclination.

Interaction Variation Maps of M i
The interaction maps of M i in Figure 11 show that M i changes in a similar manner to C i , particularly in Figure 11b. The major difference is in the stacked case in Figure 11a, where the highest M I is found when fractures are partially overlapped rather than fully stacked. When fractures are partially overlapped, the most significant shielding occurs at the center of the overlapped tips, resulting in a high value for M I . When fully overlapped, the shielding is instead uniform around the whole tip but lower, resulting in an intermediate value for M I . By measuring the area instead, C I considers the interaction magnitude to be higher in the fully stacked case and intermediate in the overlapped case, as observed in Figure 7. Figure 12 shows interaction maps that plot I , II , and III for the same geometries as in Figure 10. In the coplanar geometry in Figure 12a I demonstrates the transition between stress amplification and shielding very well. Coplanar orientations result in amplification where I = 1, and stacked orientations result in shielding where I = −1. The transition zone is well resolved even at the coarse resolution of the grid. I is presented without altering values where C I is below the threshold of interaction because the few locations where it is below 0.005 do not cause the distribution of I to be less uniform. The extent of the shielded zone changes significantly in the inclined case in Figure 12b, demonstrating the impact of fracture orientation on interaction behavior. When the interacting fracture is inclined at 45 ∘ , the extent of the shielded zone is much smaller because the region of stress amplification at the fracture tip is closest to the other fracture. At −45 ∘ , the shielded zone is much larger, because the surface of the fractures face one another in this orientation, combining two stress shielding regions. The coplanar nonoffset case of I in Figure 12a can be seen as an intermediate between these two states.

Interaction Variation Maps of i
The parameter II has the most complex distribution of the three modes. In Figure 12a, II is amplified, but a zone of shielding is found on the line x = 1 where the fractures are partially overlapped, contrasting the pure amplification when the fractures are fully overlapped. A second small zone of shielding is found when the fractures are very close. These shielding zones also arise in Figure 12b. The complex distribution of II 10.1002/2017JB014234 Figure 12. Stress intensity factor variation maps of I , II , and III for two geometries: (a) coplanar and (b) red inclined by 45 ∘ , in surface and contour format. I = 1 denotes that the SIFs around the fracture are amplified all the way around the tip and I = −1 denotes that they are shielded. Intermediate values denote that the fracture is partially amplified on one side and shielded on the other. Black dots denote sampled points, and crosses denote points where C i is less than a threshold amount (0.005). At these points, i is set to 0 because i produces scattered values when the SIFs in the isolated and interacting cases are very close. Both fractures have a radius of 1 m and 80 tip nodes.
reflects the significant changes to fracture propagation angles when fractures are en echelon, because the fracture growth angle depends on a combination of all three modes. It also corresponds to the region of lower magnitude C II in Figure 10a.
In both geometries there is generally very little variation in III . Typically, III is around 0 when the fractures are close (x and z < 2), and elsewhere the magnitude of C III is below the threshold (0.005). K III is therefore always equally shielded and amplified if it is altered due to interaction. The threshold of C i > 0.005 is mostly effective at removing misleading values where there is no interaction in modes II and III.

Comparison Between Tension and Compression
The SIF interaction results in sections 4 and 5 show results for fractures under tension. When fractures are compressed, friction between the two sides of the fracture perturbs the SIFs (Liu & Borja, 2008). This additional force is important for geological fractures, where compressive stress regimes are common and the surfaces of fractures will have significant roughness or brecciated material if they are open. In order to incorporate fracture surface friction, the gap-based augmented Lagrangian contact resolution method proposed and validated by Nejati et al. (2016) is incorporated into the simulation. This approach is specific to accurately resolving the contact and frictional forces between fracture surfaces, using an iterative approach. The method accurately calculates tractions and apertures between the two surfaces that are in contact. Initially, fractures have zero aperture, and the stick or slip condition at each fracture surface node is an emerging property of the simulation. During compression the fractures close and may either stick or slip, rendering the mode I SIFs negative until contact is resolved.
The values of SIFs depend on the friction coefficient assumed for the fracture surface. For fracture stick, SIFs are reduced to 0 in most cases, as the fracture does not display in-plane deformation. For slip cases, the deformation of the fracture leads to nonzero stress intensity factors. Figure 13 illustrates the behavior of two interacting fractures under tension and under compression. In the compressive case, the friction coefficient is assumed to be zero, in order to highlight the slip mode and accentuate interaction. Specific values of C, M i , and i for cases are as follows: 1 SIFs are plotted for the horizontal fracture. While in tension K I values dominate under compression K I values are reduced to 0, and K II and K III values become dominating. Interaction can be summarized using the same quantities as before. For tension, values of C, M i , and are consistent with low interaction. In contrast, under compression, values of C II and C III highlight the effect of the vicinal sliding fracture. In the isolated compression case, SIF values of the isolated fracture case k iso i , are close to 0. Therefore, the normalized area between SIF curves is much smaller, and the resulting C values are larger. M i captures accurately the interaction effect 10.1002/2017JB014234 between the two fractures, even if the isolated case has very low SIFs. For both tension and compression, M i yields similar values of interaction between the fractures, with the exception of mode I which is noninteractive in the compression case. It follows that the same variational maps demonstrated in previous sections can be employed to analyze the present case.
6. Discussion 6.1. Geological Context Fracture interaction maps provide a useful method for analyzing how a system of fractures will behave. Since the geometry of many fracture systems is well constrained through field and geophysical data, the state of interaction can be explored for a variety of boundary conditions. The interaction measures C and M i provide a quantitative way to establish whether fractures are interacting, by comparing the values to a threshold, as shown by their change with separation in Figure 7. For the fractures that are considered to be interacting, maps of i show whether that interaction will amplify or shield potential growth. This workflow would provide a quantitative assessment for how a network of joints or dykes will behave, for example, during coalescence into a fault or ahead of the tip of a fault with a damage zone. In particular, the amount of curvature expected when en echelon fracture pairs grow can also be analyzed using the interaction maps. Fractures growing when centered in regions where C II, III and M II, III are high will deflect significantly toward the other fracture. These regions arise when the fractures are overlapped, as found in experiments, for example, Thomas and Pollard (1993).
While the methods proposed do not vary with the scale of the fractures being studied, it should be noted that SIFs require that the fractures have distinct tips. Large faults have significant damage zones ahead of their tips, reducing how readily they can be compared to an idealized flaw. Estimations of interaction measures for faults may have to take into account the significant shear strength of faults resulting from brecciated rock inside the fault. The traction-free surface of the Earth also provides an additional stress perturbation if the scale of the fractures is comparable to their depth (Roering et al., 1997). SIF perturbations also occur as a result of heterogeneity in the medium (Yu et al., 2009). Heterogeneity is a key aspect of geological materials which can be easily incorporated in mechanical models by variation of material properties such as the Young's modulus or Poisson's ratio. This effect could be quantified using interaction measures as more general tools to study the way SIFs are perturbed. Care must be taken in order to differentiate between perturbations induced by heterogeneity, and perturbations induced by fracture interaction, because K int i will include both of these effects. This can be accomplished by calculating K int i in a volume containing all the fractures and material property variations, then calculating two versions of K iso i per fracture-one with the material property variations and one without. Then, separate values of C, M i , and i can be found, each one describing the SIF perturbation arising from either heterogeneity and interaction. Such an analysis would facilitate understanding of the competing effects of interaction and heterogeneity on growth.

Computational Considerations
The number of simulations needed to acquire a full data set with n fixed fractures in a volume is n + 1, because the system must be simulated with all fractures present, providing values for K int I-III , then once for each fracture, where that fracture is the only one in the volume, providing values for K iso I-III . Only one mesh is required for the domain, as the same mesh can be used for all calculations. The mesh does not need to be particularly refined, as the interaction measures have been shown to remain accurate for coarse meshes (Figure 8). Fracture interaction maps are relatively expensive to generate as each location of the interacting fracture requires a new mesh around the moving fracture. The method is also not restricted to using the disk-shaped domain integral method for calculating SIFs-any 3-D SIF computation method, such as displacement correlation (Branco et al., 2015) or displacement discontinuity (Thomas & Pollard, 1993), would also be effective. Variation maps can be generated using any numerical method that generates SIFs and can be also adapted to other fields.
Three-dimensional growth simulations are an additional method with potential for exploring interaction. SIF variation maps are more convenient for several reasons. SIF perturbations are reasonably computationally cheap to calculate compared to full fracture growth simulations, which may require remeshing at every step to track the propagation of the fracture. Growth simulations rely on many nonphysical parameters arising as a result of the particular modeling technique, such as the extension increment per step, which has a large effect on the final geometry of the fracture. Additionally, the geometry from each possible initial separation is difficult to summarize in a 2-D plot, unlike interaction measures which provide one numerical value per location.

Recommended Use of Interaction Measures
The interaction measures and maps detailed in this work can be used to better understand the behavior of fractures. For an outcrop or subsurface system where fracture geometries are known, examining the SIFs and interaction measures on each fracture provides a straightforward method for summarizing their behavior. Interaction maps are more useful when the locations of fractures are either variable or unknown, such as when fractures are blind, or across large outcrops. Maps shown in this work demonstrate their effectiveness at highlighting geometries which result in enhanced, reduced, or nonplanar growth. They can also be created for any pair of fractures with different geometries and orientations. Therefore, rather than referring to a catalog of maps, they should be created to fit the relative fracture orientations which are frequently encountered in a specific geological setting, following the steps outlined in section 5.
It should be noted that C and M do not directly determine the increased likelihood of propagation without additional context, such as material toughness. This can instead be approached by using a critical stress intensity factor for the material, which, if reached by K I , results in fracture propagation (Atkinson, 1984). This value can be used to scale the results of the output or in the place of K iso i to show when fractures are close enough to cause propagation. However, those results will be dependent on the properties of the medium, unlike the maps presented here, which are material and scale invariant. 6.4. Contrasting C and M i C and M i both characterize the magnitude of interaction experienced by a fracture. C yields consistent quantifications of interaction across different levels of mesh refinement. Because it is derived from comparing the area beneath two curves, its magnitude is independent of the magnitude of stress. This makes C the most convenient measure for assessing whether interaction occurs using a threshold. The benefits of using M i are that the calculation is simple; that is, it does not require numerical integration and could be found graphically, and it also provides a location on the fracture at which the SIF perturbation is highest. Either method provides similar results for how interaction changes around a fracture. In particular, the interaction map for C I in Figure 10 is almost identical to M I in Figure 11.
A benefit of using C will be when a system contains multiple close fractures. This would create a complex distribution of SIFs around the tip with multiple peaks and troughs, which would be fully captured by C, but only one peak would be selected for M i . Looking at the modes individually is also important because, in overlapped geometries, there may be multiple peaks and troughs in all three modes. Therefore, in most cases C is the recommended measure.

Conclusions
This work provides several contributions toward analyzing and understanding fracture interaction processes through study of the stress intensity factor (SIF) perturbations. These perturbations occur when two or more fractures are close enough to alter the state of stress at one another's tips. Using a 3-D finite element-based fracture mechanics model simulating pairs of fractures, these perturbations have been demonstrated and quantified. Three measures of interaction have been proposed to quantify how the SIFs are altered-circumferential SIF perturbation (C), maximum SIF perturbation (M i ), and amplification to shielding ratio ( ). C and M i can be used to quantify the magnitude of the interaction, and can be used to describe whether the stress at the tip of a fracture is mostly amplified or shielded. These measures have been used to study the interaction occurring around pairs of fractures of equal size in several orientations. This is accomplished by creating interaction maps, which plot measures of interaction in space according to the location of an interacting fracture. The analysis outlined in this work, and the interaction measures that are proposed, are applicable to any stress state or fracture population. This provides a new and powerful way to understand how fractures affect one another. tivity and the Environment (RATE) program, and European Commission for partially funding this work through the TRUST Collaborative Project, 309067. Mark Kachanov, David Pollard, Editor Yenuha Ben-Zion, and one anonymous reviewer are thanked for their constructive comments, which led to significant improvements of the manuscript. The authors also thank Véronique Lazarus for providing data for Figure 3 and Rebecca Bell and Oliver Duffy for their valuable discussions in the early stages of this work. Data presented in the graphs and interaction maps are publicly accessible through the British Geological Survey National Geoscience Data Centre (Thomas et al., 2017).