Numerical investigation of the influence of discontinuity orientations on fault permeability evolution and slip displacement

A pre-existing plane of weakness along the fault is comprised of a particular pattern of joints dipping at different orientations. The fault stress state, partially defined by the orientation of fault, determines the potential of slip failure and hence the evolution of fault permeability. Here the influence of fault orientation on permeability evolution was investigated by direct fluid injection inside fault with three different sets of fault orientations (45°, 60° and 110°), through the coupled hydromechanical (H-M) model TOUGHREACT-FLAC3D. The influence of joints pattern on slip tendency and magnitude of potential induced seismicity was also evaluated by comparing the resulted slip distance and timing. The simulation results revealed that decreasing the dip angle of the fault increases the corresponding slip tendency in the normal fault circumstance. Also, with changing joints dip angle associated with the fault, the tendency of the fault slip changes concurrently with the permeability evolution in a noticeable manner. Permeability enhancement after the onset of fault slip was observed with the three sets of fault angles, while the condition of 60° dipping angle resulted in highest enhancement. Joints pattern with a dip angle of 145° (very high dip) and 30° (very low dip) did not trigger a shear slip with seismic permeability enhancement. However, high dip and intermediate dip angles (135°, 50° and 70°) yielded high permeability in varying orders of magnitude. The large stress excitation and increasing permeability during shear deformation was noticeably high in intermediate joint dip angles but decreases as the angle increases. The magnitude of injection-induced permeability enhancement is largely influenced by the fault and joint spatial orientations. With a slight change in the joint direction, there is an increasing possibility for fault to approach a different critical state of failure. Stress elevation at the point of failure is controlled by the orientations of fault/joint planes with respect to the direction of maximum principal stress. The magnitude of injection-induced permeability enhancement is largely influenced by the fault and joint spatial orientations. With a slight change in the joint direction, there is an increasing possibility for fault to approach a different critical state of failure. Stress elevation at the point of failure is controlled by the orientations of fault/joint planes with respect to the direction of maximum principal stress.


Introduction
The quest for enhanced recovery from tight reservoirs requires a detailed study of several factors such as reservoir quality, natural fracture networks, orientation of the fractures, geomechanical properties of the matrix rock and fractures. Faults and fractures are the main targets in field development plan that enable production in naturally fractured reservoirs or induced hydraulic fractures in tight reservoirs, which practically makes fault permeability evolution study a crucial investigation in production optimization (Nelson 1985). However, this scientific study becomes more valid when variations in fault orientation and the direction of associated joints are considered in the simulation process, as this factor would strongly influence permeability anisotropy in fractured reservoirs (Watkins et al. 2018). The rapid increase in energy production has been enabled by the means of new technological advancement, such as multistage hydraulic-fracture stimulation (Rutqvist et al. 2015). Nevertheless, there are several concerns relating to the adoption of these new technologies in terms of variabilities in fault properties yielding a wide range of results (for instance, the permeability evolution as injection conditions change). Additionally, studies need to ascertain whether the injection has potential for fault reactivation, and in what magnitude is the accompanying seismicity (Davies et al. 2013;Rutqvist et al. 2013Rutqvist et al. , 2015. From reports (e.g. Shrivastava and Lawatia 2011;Xue et al. 2014;Wang et al. 2016;Feng et al. 2018;Eyinla et al. 2021), tight reservoirs have specific variations from the conventional reservoirs because of certain factors which include high heterogeneity, their deeper depth of burial, diagenetic properties, low porosities, very low permeability, poorly developed fracture system and abnormal pressure with. While rocks contain different forms of discontinuities which play an inevitable role in the overall mechanical and elastoplastic behaviour, the most significant types of discontinuities in rocks include faults, fractures, weak planes/joints, shear zones, planes of foliation, bedding planes and planes of cleavage (Eshiet and Sheng 2017). Their properties are complex, and several investigations have been carried out to assess some of their behavioural characteristics in the matrix (Brown 1987;Fairhurst 2013;Eshiet and Sheng 2017;Ghosh et al. 2018). Discontinuities influence and alter the total behaviour of rocks under in-situ and laboratory conditions. These effects are however dependent on the properties of the fault and joints, the geometry and quantity, which are also related to the locations in which they are situated in the medium. However, where two or more joints are present, the effects become more prominent, therefore, how the joints affect the behaviour of the rock is often attributed to their lower strength in comparison with the host rock (matrix) and their large-scale anisotropic properties (Eshiet and Sheng 2017).
Hydraulic fracturing via injection has become a standard technique for improving the permeability of tight reservoir in oil and gas development. However, discontinuities developed from such a process often disturb and divert hydraulic fracture propagation and path, therefore, distorting fracture fluid flow and proppant transport (Watkins et al. 2018). Consequently, the prediction of fracture/joint behaviour at varying geometry becomes an important study. Although the initiation of cracks along joint planes is mostly induced by shear failure; however, the resulting fracture reactivation is predominantly attributed to the tensile failure of the rock material (Eshiet and Sheng 2017). An essential factor which also determines the quality and yield of tight reservoirs is the distribution of fracture-controlled permeability resulting from fluid injection, which could be attributed to several factors such as the pattern of joints. Nonetheless, one way of characterizing these fractures is by the application of numerical forward modeling (Parker 2013) and laboratory experiments (Asahina et al. 2018;Feng et al. 2018). The combination of both would give a more detailed study and allows adequate scientific correlation. Field developments now involve the adoption of improved static reservoir characterization method, which incorporates essentially the geomechanical properties of the field and the initial stress distribution on the reservoir, including the numerical reservoir modeling which examines the dynamic evolution of stress state during fault injection (Turner et al. 2017).
Studies have shown that the spatial distribution and orientation of fault may or may not be the same as the associated joints (Cappa and Rutqvist 2011), but the existence of these discontinuities in different forms have an intrinsic influence on the numerical simulation for thermo-hydro-mechanical interactions. On a wider scale, it determines the expected hydraulic fractures emerging, and the overall recovery process (Men et al. 2018). Many numerical and analytical solutions to both hydraulic fracturing problem and Enhanced Geothermal Systems have been proposed, and each of these has improved the understanding of the thermo-hydro-mechanical response of fault under injection, especially when it considers the influence of fracture geometry (Adachi et al. 2007;Rutqvist et al. 2013Rutqvist et al. , 2015Yang and Zoback 2014;Jacquey et al. 2015;Gan and Elsworth 2016a;Feng et al. 2018). Additionally, many studies have reported field experiments and laboratory studies of hydraulic fracture behaviour in both large and small scales (Blair et al. 1989;Legarth et al. 2005;Roberts 2005;Casas et al. 2006;Athavale and Miskimins 2008;Roberts and Abdel-Fattah 2009;Liu and Manga 2009;Chuprakov et al. 2010;Elkhoury et al. 2011;Faoro et al. 2012;Candela et al. 2014). From these reports, permeability enhancement mechanism is characterized by dynamic stresses induced by fluid pressure.
The study conducted by Gan and Elsworth (2016a) explored a diverse stimulation scheme to determine the impact of stimulation direction relative to the fracture orientations on the magnitude and extent of thermal recovery rates, for a proposed reservoir with a defined pre-existing fracture network. However, to create a discrete fracture network in a simulation, several factors are often considered, which include, fracture location, the orientation, length, and fracture aperture (Gan and Elsworth 2016b). Cappa and Rutqvist (2011) reported that the most important factors for the initiation of fault slip are shear stress development and fluid pressure, as they are known to enhance fault rupture by overcoming the fault frictional resistance. Consequently, for any numerical study, slip-tendency analysis provides a technique which allows scientific evaluation of stress states and how it can be related to expected seismic or aseismic activity. Resistance to frictional sliding has been identified as a factor which is responsible for the slip behaviour of faults during fluid injection (Jacquey et al. 2015), and it is influenced by the properties of the associated joints (Zhang et al. 2018).
Generally, during unloading of fault (reduction in the normal stress), hysteresis in the fracture permeability is often observed. This implies that fault permeability is usually noticeably different during unloading from the loading phase even when they are stressed equally (Gutierrez et al. 2000). The hysteretic property of the fault permeability under normal loading simply connotes a manifestation of the general irreversibility of rock deformation (Lavrov 2017). However, concerning the simple relationship between fractures, stress and strain during injection, contributions from many studies have presented theories, numerical approach and analysis of injection-induced fault reactivation and seismic slip (e.g. Adachi et al. 2007;Rutqvist et al. 2013Rutqvist et al. , 2015Yang and Zoback 2014;Elsworth 2014, 2016b;Jacquey et al. 2015). These and many more reports have enhanced the understanding of the fault reactivation processes, especially in enhanced geothermal systems. However, because of the role played by fault orientations in the stress field analysis, variations in fault permeability with respect to the changing fracture/joint orientations can now be explained through numerical simulations. An example of these studies was earlier presented by Jacquey et al. (2015), demonstrating how the angle of fault influences the initial slip tendency and dynamic permeability evolution. It is crucial to understand how permeability changes during shear deformation, as this would afford an understanding of the dynamic hydromechanical processes during injection and mechanisms influencing the occurrence of earthquakes both at shallow and in deep crustal levels (Tanikawa et al. 2010).
Thus, this study explores the roles of fault geometry and associated joints, using data from tight shale reservoir from Akas field of the Niger Delta Basin as a case study. Studies from this field have been discussed in previous reports (Eyinla and Oladunjoye 2019;Eyinla et al. 2020Eyinla et al. , 2021. The reservoir model is presented as a finite medium with a hydraulically induced normal fault, and the overall mechanical behaviour of fault is represented by a set of solid elements with ubiquitous joints which are oriented as weak planes in the fault zone as described by Cappa and Rutqvist (2011). We investigate how their relationship could modify the poroelastic response of the fault under undrained simulation conditions. In overall, the study seeks to understand the response of fault slip behaviour during injection as orientational properties change for the purpose of production optimization, reservoir management and prediction of seismic event.
2 Theory of study and methodology

Stress theory
The effective normal stress and the shear stress acting on a fault ( Fig. 1a) are estimated from the generated data after simulation, using the approach modified after Gan and Elsworth (2014) with the parameters in Fig. 1a thus: where r n represents the effective normal stress, P is the pore fluid pressure, r 1 is the maximum principal stress, r 3 is the minimum principal stress and h represents the angle between the fault plane and the maximum principal stress, r 1 direction. Therefore, Coulomb stress ratio, g, which is defined as the ratio of shear stress to effective normal stress (Fig. 1b) is represented as: Notably, for any surface, the slip tendency can be described as the ratio of acting shear stress s ð Þ to the effective normal stress r n ð Þ. However, a slip is likely to occur on the surface when this ratio is greater than or equal to the frictional resistance to sliding (Jaeger et al. 2009). Also, the static friction coefficient, l s ; has been defined by Biot (1941) and Byerlee (1978) as: where / is the friction angle. And for this study, the internal fault friction angle used is 28°, consequently, the coefficient of friction, l s = 0.53. Therefore, for a slip to occur, the maximum coulomb stress ratio value (called the peak friction) must be greater than or equal to 0.53 ðg ¼ s=r n ! l s Þ. Thus, the reactivation of a pre-existing fault is likely to occur during fluid injections, depending on the maximum sustainable pressure limit and principal stress resolution (Kim and Hosseini 2014).

Fault permeability and aperture evolution
In this model, a demonstration of the sensitivity of fault permeability to hydromechanical interactions, normal stress change as well as volumetric strain is presented. The behavior of the fault would undoubtedly influence a change in the fault normal displacement. Consequently, sudden increase in the fault permeability would result at the onset of shear slip (Eyinla et al. 2020). For a fractured medium, models for permeability change as governed by the input variables involves the growth in the fracture aperture, which may be defined by an exponential function of applied effective stress r and the nonlinear fracture stiffness a (Rutqvist et al. 2002; Elsworth 2014) given as: where b is the current hydraulic aperture due to current effective normal stress, r 1 n , br is the residual aperture, b max is the maximum aperture without mechanical stress effect, r 1 0 is the effective stress at which zero deformation occurs (usually 0), a is the non-linear fracture stiffness. Generally, the permeability of the fault damage zone is better presented when the pattern of the fracture is well defined, which however requires assigning the right aperture and fracture spacing. As a result, the initial permeability of the fault zone is higher than the permeability of the matrix (Table 1), and the transmissivity of fluid pressure within the fault zone is related to the hydraulic aperture of the fracture (Norbeck and Horne 2015). However, to satisfy the requirement for representing a coupled non-linear elastic behaviour of fault, the permeability of fractures in the fault zone has been modelled using existing approach (Warren and Root 1963) as it connects the fracture aperture and fracture spacing through the relation: where k is the fault permeability (m 2 ), b is the fracture aperture (m), and s is the fracture spacing (m).

Mechanism of shear failure and seismic slip
The fundamental mechanism of fault reactivation is expressed when the shear stress exceeds the shear strength of the fault (Cappa and Rutqvist 2011). Slip tendency is therefore regarded as the likelihood of a surface to slip during injection, a mechanism which is highly dependent on the frictional resistance of the fault and measured by the ratio of shear to normal stress acting on the fault plane ( Fig. 1b). Analysis of stress distribution on the fault plane relative to the shear strength provides an understanding of its potential to cause a slip, and thus, makes it possible to evaluate general exploration risks, including seismicrisk, fault-rupture risk assessment and earthquake forecasting (Morris et al. 1996). The commonly used relationship describing fault slip in the failure analysis Thermal expansion coefficient of solid (K À1 )] 12 9 10 À6 12 9 10 À6 12 9 10 À6 And according to the Terzaghi (1923), the effective stress is expressed as: where s is the critical shear stress for slip occurrence, c is the cohesion, l s is the static friction coefficient, r 0 n is the effective normal stress, r n is the total normal stress and P is the fluid pressure. An approach to distinctly determine the fault stability and if a failure has occurred on a fault plane dipping in the vertical direction has been presented by Rutqvist and Oldenburg (2007) and Jaeger et al. (2009). This is examined by comparing the ratio of the maximum and minimum principal effective stresses with the frictional resistance as: Here, q represents the effective stress limiting ratio according to Biot effective stress theory (Biot 1941), and is defined as, where l s ¼ 0:53, (from Eq. 6), the value of q ¼ 2:76.
In the initial setting of the model, a is the Biot coefficient sets at 1.0, r 1 = 45.5 MPa, r 3 = 27.3 MPa and P = 13.8 MPa, therefore, r 0 1 =r 0 3 ¼ 2.3. Consequently, the initial status of the fault is stable since 2.3 \ q = 2.76.

Estimating seismic magnitude
Magnitude of seismic events in the fault zone as a result of fluid injection can be quantified using seismological theories and adopting the approach of Cappa and Rutqvist (2011) and Mazzoldi et al. (2012). This is based on a moment magnitude scale which describes the strength of the seismic event according to the energy released by the seismic slip in the fault plane (Kanamori and Abe 1979). The first step is to quantify the size of the seismic event after simulation for the ruptured surface of the fault zone. This attribute has been described as the seismic moment, M o ; defined by Kanamori and Brodsky (2001) as: From Aki (1967), this expression can also be rewritten as: Here M o is the seismic moment Nm ð Þ, l is the shear modulus Pa ð Þ, A is the rupture area (m 2 Þ, D c is the mean slip m ð Þ, L is the fault length m ð Þ, and W is the fault rupture width m ð Þ. Estimating the moment magnitude M ð Þ of the seismic event involves the adoption of an equation which relates seismic moment, as given by Kanamori and Anderson (1975) thus: This M À M o relationship can also be expressed (Kanamori and Abe 1979;Purcaru and Berckhemer 1982) as: 3 Model analysis The simulation for this study involved the use of the coupled hydro-mechanical simulator TOUGH-REACT-FLAC3D developed by Taron et al. (2009), which links the TOUGHREACT multiphase flow with the FLAC3D geomechanical simulator (Itasca 2009). The elastoplastic behaviour of the fault in FLAC3D which occurs as a ubiquitous fractured media impressively represent an anisotropic mechanical behaviour.
A coupled hydromechanical fault model can be developed within the framework of TOUGH-FLAC by utilizing existing capabilities within TOUGH2 and FLAC3D codes, and by developing specially designed coupling modules for faults. The fault model in TOUGH-FLAC could be discretized in terms of the mechanical behaviour of faults and fault zones represented in FLAC3D by either special zero-thickness mechanical interfaces (Fig. 2a), by an equivalent continuum representation using solid elements (Fig. 2b), or by a combination of solid elements and ubiquitous-joints oriented as weak planes (Fig. 2c).
One merit of using solid elements as ubiquitous fractured media is the ability to model faults to account for parallel and cross fault heterogeneity. The mechanical behaviour of the fault as a ubiquitous fractured media accounts for the presence of an orientation of weakness (weak plane) in Mohr-Coulomb model (Cappa and Rutqvist 2011). The Mohr-Coulomb envelope has a tension cut off which serves as the criterion for failure on the fault weak planes (Hacker 1997). Using this procedure, it is possible to model the plastic flow behaviour for both the weak planes and rock matrix in the vicinity of the fault zone.
The fault is designed to contain a low permeability fault core which is flanked by slightly higher permeability damage zones. The elastoplastic constitutive model available in the FLAC3D was adopted in the study. The matrix was assigned Mohr-Coulomb model, while the constitutive model for the major fault was defined as strain-hardening/softening ubiquitous (subi) joint model. This choice is adequate to model the planes of weakness introduced by the fracture zones. The fault is set to be critically stressed, dipping towards the direction of the maximum principal stress (Fig. 3a). The fault architecture is designed with finer mesh than the other part of the reservoir. That is, the mesh size in fault and the matrix to the left and right of the fault zone contains uniform and smaller sizes than those in the other upper and lower regions of the matrix (Fig. 3b-d) to ensure accurate and efficient simulation of the zone of interest. The model was constructed using structured block grids in FLAC. Table 1 shows the assumed material properties for the fault zone and host rock, derived from laboratory measurements and previously published data from the study area (Eyinla and Oladunjoye 2019;Chukwu 2017;Nwozor et al. 2017;Emudianughe and Ogagarue 2018;Ichenwo and Olatunji 2018;Ogunsakin et al. 2019;Eyinla et al. 2020Eyinla et al. , 2021. The constitutive mechanical properties for the model were derived from the work of Gan and Lei (2020). The study presents a simple reservoir model geometry with spatial dimension of 600 m 9 15 m 9 600 m (x, y, z), including a single hydraulically induced normal fault integrated in the model at interval 150 to 450 m (z-direction). For each model, three distinct fault orientations were examined, dipping at NE 45°, NE 60°and NW 20°denoted as angle 110° (Fig. 3b-d).
The fault is represented as solid elements with ubiquitous joints existing as weak planes with specific orientations with respect to the strike of the fault plane. Figure 3a is a schematic representation of 45°fault angle, showing the assumed associated joints and their orientations. In this study, the specific joint direction assigned to each fault orientation are 145°, 135°, 70°, 50°and 30° (Table 2). The model geometry constructed consists of 1320 elements: 1080 elements for the matrix, and 240 elements for the fault zone. These zones were populated with material properties in Table 1. This include hydraulic properties (e.g., porosity, permeability), elastic moduli (e.g., shear modulus, bulk modulus), material strength properties (e.g., tensile strength, friction angle, cohesion), amongst others. The elastoplastic constitutive model available in the FLAC3D was adopted in the study. The matrix was assigned Mohr-Coulomb model, while the constitutive model for the major fault was defined as strainhardening/softening ubiquitous (subi) joint model. This choice is adequate to model the planes of weakness introduced by the fracture zones. The fault architecture is designed with finer mesh than the other part of the reservoir. That is, the mesh size in fault and the matrix to the left and right of the fault zone contains uniform and smaller sizes than those in the other upper and lower regions of the matrix (Fig. 3bd) to ensure accurate and efficient simulation of the zone of interest. The model was constructed using structured block grids in FLAC.
It is expected that the fault core would serve as a barrier preventing fluid from penetrating across it, whereas the damage zones create permeable channels through which fluids are transmitted in parallel direction to the strike of the fault. Model boundaries are set at no flow boundaries, and the y-direction is set at roller condition with no normal displacement as constant stresses are applied (Gan and Lei 2020). The initial permeability of the fault core and damage zone are sets at 10 -15 m 2 and 10 -14 m 2 respectively, while the initial matrix permeability is set at 10 -16 m 2 . An interesting component of this model is the ability of the fault permeability to evolve through time. Therefore, at the onset of slip, it is expected that the fault zone permeability would change, and there would also be a drastic reduction in the shear stress.
The friction angle of the fault joint is 28°while the dilation angle is 5°. The injection well is centrally located along the fault mid-point and extends laterally towards the y axis (15 m), with a constant injection rate at 0.3 kg/s. The injection condition for this study is isothermal. The minimum principal stress acting on x-direction was 27.3 MPa and the maximum principal stress was set at 45.5 MPa (z-direction), they are illustrated in Fig. 3a. The permeability evolution model in FLAC3D adopted in this work has been verified in earlier studies by White et al. (2016). The behaviour of the model promotes permeability enhancement at the onset of fault slip as a result of the dilation angle of the fault.

Simulation scenarios
In the first scenario (Scenario A), there are three basic fault configurations (45°, 60°and 110°) with the same direction of joint which was assumed to be 135°as shown in Table 2. Here we examined the sensitivity of each of the fault angles to hydromechanical processes induced by injecting fluids at a constant flow rate of 0.3 kg/s. We then further explored the effects of changing joint directions on the HM response of the three fault angles under the 5 chosen joint directions (30°, 50°, 70°, 135°and 145°) with the same injection b Fig. 3  rate at 0.3 kg/s. This investigation is to fully examine and certify the effect of each joint direction on the hydromechanical behaviour of the individual fault orientations. This is Scenario B, and the variations in fault/joint orientations was examined for permeability evolution as well as the magnitude of injectioninduced seismicity under the same matrix-fault properties and injection conditions. For the slip displacement measurement, points along the fault surface were recorded from A to G. Point A being the top of the fault and point G is the base (Fig. 3e).

Results and discussion
An investigation into the effects of the frictional property and compaction dominating behaviour of the fault in this study has shown that HM interactions in the fault zone change with different fault angles. The results of the three fault angles 45°, 60°and 110°for Scenario A (as shown in Table 2) revealed that angle 60°has the highest permeability enhancement (Fig. 4). The timing of fault failure with changing fault angle is earliest in 45°and most delayed in 110° (  Fig. 4a-e). Effective stress and shear stress magnitudes are greater at lower fault orientations as 45°[ 60°[ 110° (Fig. 4d-e), and stress drop was highest at fault angle of 45°. However, the simulation result after changing the joint direction for each fault orientation showed how the joint greatly influence the fault stability and how the impact on permeability evolution for each model was revealed (Figs. 5, 6, 7). Also, the variation in shear failure tendency and eventual slip timing is different, as the slip potential is highly dependent on the joint direction with respect to the existing fault plane. This implies that with a change in the joint direction, there is an increasing possibility to have a different levels of fault criticality to failure. Additionally, Fig. 8 shows the correlation of each fault/joint model as shear slip is initiated when the shear strength limit (coefficient of friction,l s ¼ 0.53) is exceeded, and the principal effective stress ratio ðr 1 0 =r 3 0 Þ is greater than the effective stress limiting ratio (q ¼ 2:76Þ. A slip will probably not occur even when the assumption about fault angle relationship with the direction of maximum stress is fulfilled if the angle between the joint and the principal stress is very large, just as is the case of 145°and 30°j oint orientations, the fault tends to be relatively stable. The slip profile for the three fault angles under various joint configurations illustrate the magnitude of slip events much better in Figs. 9, 10, 11 and 12. The slip curves are dissimilar for each fault/joint model and the magnitude of fault displacement depends on the joint direction in relation to the fault plane and the direction of maximum principal stress.

Effect of fault orientation
Variation in the orientation of the fault dip created significant influence on the loading response, permeability evolution (Fig. 4a) and the slip tendency of the fault. It was observed that increasing the dip angle of the fault delays the slip timing (Fig. 4b), however, the three orientation angles displayed tendency of producing fracture slip under same hydromechanical properties and injection conditions adopted for the simulation. The onset of slip timing increases as the fault angle decreases. In this way, fault configuration with angle 45°resulted in an earlier onset of slip after 15.7 days of injection, whereas slip occurred only after 16.7 days of injection in the case of 60°, and 19 days for fault angle 110°. This implies that the slip tendency of a pre-existing fault is controlled not only by the strength of fault (Eshiet and Sheng 2017;Eyinla et al. 2020), but also by the orientation of the fault. Also, the magnitude of permeability enhancement after the onset of slip is significantly influenced by the fault orientations, as fault permeability was enhanced by one order of magnitude in angle 110°(from 10 -15 to 10 -14 ), two orders of magnitude in the case of angle 45°(from 10 -15 to 10 -13 ), and over four orders of magnitude in angle 60°(from 10 -15 to 10 -10 ). From the work of Akande et al. (2021) that involved the same permeability evolution model, 60°fault dip angle yielded the most interesting result for their study after comparing results from other fault angles. This is obviously the same with our study, where fault dip angle 60°resulted in the highest permeability when compared with the other two fault angles.
The tendency of permeability enhancement as injection progresses arguably results from fault opening during sliding, mostly reflecting the roughness of the fracture walls and the effect of the associated microcrack dilatation occurring during the early stages of the fault failure (Guglielmi et al. 2015). However, the peak behaviour of each fault angle is well represented on the coulomb stress plot in Fig. 4b. The plot illustrates variation in frictional resistance of the fault as the orientation changes, despite the controlling effect of the peak friction angle and the residual friction angle (Gan and Lei 2020). Previous reports by Cappa et al. (2018) and Gan and Lei (2020) showed how faults respond in terms of permeability enhancement when fault reactivation occurs. Seemingly, the Continuous permeability elevation after shear slip is evident in all the fault angles, with continuous pore pressure build-up. A steadily rising pore pressure is the usual poroelastic response of faulted reservoir to fluid injection, regardless of the injection fluid properties (Altmann et al. 2010;Schoenball et al. 2010;Vilarrasa et al. 2013;Kim and Hosseini 2014). Similarly, our results showed a continuous pressure build-up along the fault plane, meanwhile, the variation in this hydraulic response is a function of the aperture and the fault initial stability. This stability is dependent on the fault dip angle, which determines the angle between the fault plane and the direction of principal stress orientation (Streit and Hills 2004). There is a drop in pressure at the point of slip observed for fault angles 110°and 60°due to the significant void space creation through shear dilation. Conversely, there is a sudden rise in pressure just before the onset of shear failure for the case of angle 45°. It might be caused by the compaction effect of the pressure build-up in undrain conditions, with the increasing normal stress (Fig. 4d).
Clearly, there is an intricate interaction between fluid pressure, fault permeability and the corresponding deformation (Cappa et al. 2018). The result for the three fault angles showed that the fault permeability evolution is also pressure dependent, considering the increasing pressure observed as permeability enhancement evolves, and the corresponding sudden pressure excitation/relaxation as permeability changes with fault slip. This observation is similar to earlier reports by Barree et al. (2009), Cho et al. (2013, Wang and Sharma (2019), and Gan and Lei (2020). However, there is a higher stress sensitivity for low permeability fractured reservoirs, which obviously diminishes when the effective stress rises above certain level (Archer 2008). This study observes a steady drop in stress as simulation commences (Fig. 4d), yielding a gradual permeability rise before slip (Fig. 4a).
Notably, when the fault plane reached a critical stress state, the unloading process was accompanied by a released strain energy. The fault angle with more shear stress produces earlier onset of slip, larger shear drops and larger slip zone, and this finding is consistent with reports in earlier studies by Rutqvist et al. (2015) and Guglielmi et al. (2015). This revealed why fault angle 45°which results in a larger shear accumulation among the three fault orientations considered in this study also produced the earliest onset of slip, the largest shear drop, and the largest slip area. Although the slip distance in fault angle 60°is highest among the three cases (Fig. 9a), the overall slip area in angle 45°is the largest with 5 slip points (Fig. 9b), and this is a direct reflection of the magnitude of shear stress drop (Fig. 4e). This means that a fault with higher background stress would produce larger slip growth simultaneously (Cappa et al. 2018). The stress reduction during unloading is mostly due to interparticle force decrease and continuous particle contact breakage. Nonetheless, the effect of changing frictional property and compaction dominating behaviour of the fault as orientation changes has greatly influenced the variation in shear slip timing.
Furthermore, in response to loading, while permeability is significantly enhanced when stress is elevated, effective stress evolution affects both the timing of the slip and the size of the corresponding slip area. Meanwhile, stress anisotropy is increasingly developed as the angle between the fault plane and the maximum principal stress increases. Fault angle 45°b eing the highest in this case developed the highest magnitude of effective stress, followed by angle 60°. However, regardless of the effective stress and shear stress elevation and magnitudes, the resulting final permeability enhancement after shear failure is highly dependent on the shear strain and shear dilation of the elasto-plastic model as reported by Rutqvist et al. (2013).

Effect of joint orientation
The previous section assumed 135°as the joint dip angle for the three fault orientations, indicating that the results obtained did not consider the effect of the  Figure 5a-d, however, show the effect of various joint dips on the fault permeability enhancement among the three fault orientations. The fault configurations are plotted together to discern the variation of fault permeability in the same joint direction. Although the relationship of the fault plane with respect to the maximum principal stress direction has been observed to determine fault stability, our results have showed that fault slip potential is also dependent on the joint direction with respect to the existing fault plane (Figs. 5, 6, 7). Literally, a slip will probably not occur when the earlier assumption about fault-maximum stress direction relationship (12, 14, 39) is fulfilled. That is, if the angle between the joint plane and the maximum principal stress direction is large, as in the case of joint dip 145°and 30°in Fig. 5a and d. The effect of joint direction in these two cases promote very high fault compaction and increasing frictional resistance which influences the fault plane to resist shear deformation. This is because with the extremely high stability of the fault plane induced by the joint, fault slip tendency is negligible, thus, no amount of fluid pressurization is sufficient to overcome the frictional strength and the resistance to slip within the limit of our experiments. Overall, as the stability becomes more enhanced in the fault plane, fault slip tendency is lowered. A related observation was reported by Grasselli and Egger (2003) and Jacquey et al. (2015).
As this study has revealed, the most significant fault permeability enhancement is only observed when there is a shear slip. From result, the variation in slip Evolution of pore pressure acting on transition near fault centre with the three fault angles and joint dip 50°1 tendency and permeability evolution is enhanced as the joint angle changes indicating that permeability enhancement is strongly affected by the joint orientation in the model, as it affects the stress state and slip potential. Here, joint direction of angle 50°produced the highest permeability enhancement for the three fault angles (Fig. 5c), invariably, the highest slip distance in each fault/joint model. This implies that with a changing joint direction, there is an increasing possibility to have a different levels of fault criticality to failure. In the case of fault angle 110°, all the associated joint dip angles generally produced low permeability enhancement compared to the other two faults angles. The enhancement after slip lies within two orders of magnitudes, even with intermediate joint dip 50°, yielding the most significant permeability increase for fault angle 110° (  Fig. 6a). However, in fault angle 45°, joint dip 50°a lso produced the highest permeability enhancement, but it is quite more enhanced than it was in fault angle 110°, increasing from the order of 10 -15 m 2 to 10 -12 m 2 . (Fig. 6c). This is higher than the enhancement observed in fault angle 45°with joint dip 135° (  Fig. 4a), which increases from the order of 10 -15 m 2 to 10 -13 m 2 , whereas 45°joint dip 70°only increases with less than an order of magnitude after slip (from 10 -15 m 2 to 10 -14 m 2 ). Nevertheless, fault angle 60°with joint dip 135°produced the highest permeability enhancement (Fig. 4a) when compared to the others under CASE C ( Table 2, and Fig. 6b). This also vividly suggests an increasing slip distance during shear failure.
Evidently, as the joint direction changes, the mechanism of slip reactivation on preexisting faults when the Coulomb failure point is reached respond to the changing frictional resistance. The resultant effect ig. 12 Fault slip variation of fault angle a 60°with different joint dip orientations. b 45°with different joint dip orientations c 110°with different joint dip orientations of this is the observed variation in Coulomb failure point, where the variation in critical peak friction value highlights the effects of the joint dips (Fig. 7a,  c). However, the steep pore pressure surge noticed in Fig. 8b and d reflects the timing when the shear slip occurs. Interestingly, we observed a less pronounced slip event in the case of joint dip 70°as shown in Figs. 10c and 12. The joint plane in this case is closer to the vertical (maximum stress direction). However, intermediate joint orientation (50°and 135°) causes friction weakening which promotes lowering of fault strength, consequently, influencing accelerated slip (Figs. 9,11,12). In this study, it is observed that the rate at which fault reaches its peak strength is much higher with an intermediate joint dip and fault orientation (Fig. 8).
Nevertheless, since earlier studies using the same model and simulation parameters have underscore the response of the fault under different injection position and flow rates (Eyinla et al. 2020), it will be worthwhile to investigate thermal influence on the behaviour of each fault/joint model in terms of fault reactivation potential and seismicity magnitude. Moreover, for a complex fracture system, the combined effect of the different orientations might produce a significant contribution to the overall result.

Slip distance analysis
To further investigate the extent and magnitude of shear failure as fluid pressure induces frictional behaviour, this study explores the slip distance measurement to illustrate the impact of fault and joint orientations on induced seismicity during fluid injection. Notably, during injection, critical high fluid pressure causes fault opening and seismic slip when the shear strength of the fault is exceeded. However, the overall effect of the injection process is mostly felt around the injector location (Cappa et al. 2018), and the induced seismicity, which is created by the continuous fluid injection, is often activated indirectly at some distance away from the injector position via stress transfer mechanism associated with disseminating fault slip propagation (Guglielmi et al. 2015). This model showed that enhanced permeability favours the growth of aseismic slip even beyond the pressurized area. The seismic rupture propagates away from the injector to enhance permeability through shear dilation, although the weakening/reduction of fault strength is more pronounced at the pressurized zone than in the immediate vicinity. However, when monitoring near the injection point, where fluid pressure is most elevated, aseismic slip is predominantly driven by effective stress reduction.
The slip magnitude is also directly linked to the growth of permeability, the larger slip distance leads to the higher permeability enhancement in the scenario of 60°dip angle. (Figs. 4,9). In this case, the slip distance profile along the fault damage zone shows how the size of slip zone varies as the fault angle changes. Fault angle 60°delivers the highest slip distance in the range 0.025-0.045 m in all the three joint dip orientation (Figs. 9-12). This is followed by fault angle 45°with the highest slip distance being recorded in joint direction 50°as 0.16 m (Fig. 11b). Generally, fault angle 110°delivers the lowest slip distance, regardless of the associated joint orientation (Figs. 9c,10c,11c,12c). The slip magnitude in the case of fault angle 110°is highest in joint dip 50°with slip distance 0.012 m (Fig. 11c). This therefore implies that under same injection conditions, the most unstable among the three fault orientations is 60°, with the highest likelihood of a slip at most joint orientations (Figs. 9a, 10a, 11a, 12a).
There is a higher fault permeability yield when the fault dip angle has lower frictional resistance, thus, producing high slip distance. Thus, increased shear stress produces faster onset of slip, and this results in a larger slip area (Cappa et al. 2018), which was observed in the case of fault 45°with highest shear stress (Fig. 4e) and producing the largest slip area (Figs. 9b,11b,12b). The influence of joint dip angle on slip distance and permeability evolution was also evident, as angles with intermediate joint dips are observed to induce high slip displacement (Figs. 9, 11, 12) than very high dip angles (Figs. 10, 12).
It was believed that mesoscale fractures and inter and intragranular cracks are developed within fractures with high slip displacement (Tanikawa et al. 2010). Apparently, this established why fault angle 60°with joint dip 135°produced the highest permeability enhancement at the end of injection, from initial order of 10 -15 m 2 to 10 -10 m 2 compared to the other two fault angles. Thus, permeability enhancement during fault injection varies by different orders of magnitude (generally from 1 to 5), depending on the fault/joint orientation. However, the induced seismicity on the fault plane for the various fault/joint orientations have moment magnitude, M, ranging from 0 to 1.0, indicating very low aseismic event.
Here, seismicity is highest in fault angle 60°jdip 135°t han the other two fault/joint models, which suggests that the moment magnitude is dependent on the fault initial stability being proportional to the associated joint direction. Thus, the amount of slip resulting from fault failure determines the magnitude of seismicity.

Conclusions
This paper explores how the frictional equilibrium of pre-existing fault is altered in a diversified range as stress redistribution in the vicinity of the fault changes. Distinctly, the mechanism driving unloading-induced fault instability have been presented with much emphasis on the frictional resistance of the fault due to changing orientation using the Mohr-Coulomb model. This study has presented an exclusive description of how different fault orientations and associated joint directions may interact with the stress distribution during direct injection into the fault zone for the purpose of enhancing the permeability. Fault reactivation and injection-induced seismicity, in terms of timing and slip displacement magnitude, have been found to depend on different fracture configurations. Evidently, fault reactivations occurred when the accumulated shear stress acting on the fault plane exceeds its shear resistance, and the effective stress in the fault zone must respond to loading before any fault failure can occur.
The following summary resulted from the fluid injection simulation study: 1. Lower fault angles generally favor early onset of fault slips (45°[ 60°[ 110°), but the absolute effect of fault configurations is further modulated by the directions of associated joints with respect to the fault plane and the direction of maximum principal stress (r 1 ). That is, both fault and associated joints' direction have a first-order control over permeability evolution and fault reactivation during fluid injection. 2. With increasing distance of the joint plane from r 1 ; frictional strength of the joint increases, and the joint planes thus promote further fault compaction. For instance, joints with very low dip oriented at 030°and very high dip oriented at 145°r esult in the lowest and insignificant fault permeability enhancement during the stimulation in all fault scenarios. These two joint directions yielded no fault failure because the effect of their direction confer special frictional stability on the fault plane instead of inducing fault slip as in other joint directions examined (50°, 70°, and 135°). 3. The slip timing and slip tendency change as the fault angle and joint orientation change. This result is also evident in the permeability evolution. However, joint direction of angle 50°produced the earliest onset of slip in all fault angles. It also resulted in the highest permeability enhancement for the three fault angles, and invariably, the highest slip distance among all the simulation scenarios of fault/joint model. 4. Permeability enhancement during fault injection varies in order of magnitudes depending on the position of fault plane relative to joint orientation, and the incurred slip distances have positive contribution in enhancing the magnitude of joint aperture through shear dilation, thereby enhancing the permeability of fault after slip. 5. Both the effective stress and shear stresses have a steep rise when the fault approaches a critical state of failure, and the magnitude of this elevated stresses is highly sensitive to the orientation geometry of the fault and joint. Fault angle with higher effective stress during unloading is also consistent with developing higher shear stress, as seen in the case of angle 45°. Meanwhile, this similar order of effective stress and shear stress magnitude define the timing and onset of slip, where fault angle 45°yielded the largest shear and effective stress magnitude and thus results in earliest onset of slip. 6. The seismicity resulting from fluid injection depends on the slip tendency and the slip magnitude which has a direct correlation with the fault initial stability in terms of the relationship with the associated joint. The fault angle which produced the highest shear enhanced permeability produced the largest aseismic event, and from our study, fault angle 60°jdip 135°produced the most enhanced permeability, up to 10 -10 m 2 and thus, the highest slip displacement and seismicity (M = 1.0).
The implications from this study would enable a more accurate injection well placement and guidance when exploiting faults and understanding their behaviours as conduits during recovering subsurface energy so that production can be maximized. There is a direct connection between permeability changes and effective stress changes in the fault zone, which is also proportional to the magnitude of accompanied seismic event and slip zone.