Relationship Between the Orientation of Maximum Permeability and Intermediate Principal Stress in Fractured Rocks

Flow and transport properties of fractured rock masses are a function of geometrical structures across many scales. These structures result from physical processes and states and are highly anisotropic in nature. Fracture surfaces often tend to be shifted with respect to each other, which is generally a result of stress‐induced displacements. This shift controls the fracture's transmissivity through the pore space that forms from the created mismatch between the surfaces. This transmissivity is anisotropic and greater in the direction perpendicular to the displacement. A contact mechanics‐based, first‐principle numerical approach is developed to investigate the effects that this shear‐induced transmissivity anisotropy has on the overall permeability of a fractured rock mass. Deformation of the rock and contact between fracture surfaces is computed in three dimensions at two scales. At the rock mass scale, fractures are treated as planar discontinuities along which displacements and tractions are resolved. Contact between the individual rough fracture surfaces is solved for each fracture at the small scale to find the stiffness and transmissivity that result from shear‐induced dilation and elastic compression. Results show that, given isotropic fracture networks, the direction of maximum permeability of a fractured rock mass tends to be aligned with the direction of the intermediate principal stress. This reflects the fact that fractures have the most pronounced slip in the plane of the maximum and minimum principal stresses, and for individual fractures transmissivity is most pronounced in the direction perpendicular to this slip.


Introduction
The permeability of subsurface fractured rock masses is important in the context of natural aquifer dynamics, storage and retrieval of gases and fluids from the subsurface, geological waste disposal, and tunnelling and mining operations. In most cases, and in particular when the matrix rock is not highly permeable, the overall fluid conductivity of the fractured rock mass is strongly influenced by fractures.
Anisotropy is intrinsic to fractured rock permeability, for various reasons. The state of stress at any point in the Earth's crust is often anisotropic, with a significant deviatoric component (Zoback et al., 1989). Fracture networks that form under these conditions are often geometrically anisotropic and display highly directional patterns that often record the deformation history of the rock mass (Pollard & Aydin, 1988). Depending on their individual orientation, fractures will experience different magnitudes of compression and shear displacement, both of which have a strong influence on hydraulic transmissivity. In some fractured rock masses, the orientation of the global permeability tensor may be controlled more by dissolution/precipitation processes than by the remote stress field, as has been pointed out by Laubach and his collaborators (Laubach et al., 2004). The present paper does not address this scenario but rather focuses on those rock masses whose permeability is strongly controlled by the current stress state. Moreover, fractures having significant shear displacement will have locally anisotropic transmissivity, with higher transmissivity in the direction perpendicular to the shear direction than parallel to it, as has been shown theoretically (Mallikamas & Rajaram, 2005) and experimentally (Auradou et al., 2005;Gentier & Hopkins, 1997;Matsuki et al., 2006;Nemoto et al., 2009;Yeo et al., 1998). This occurs primarily because aperture correlation is less persistent in the sliding direction as compared to the normal direction (Mallikamas & Rajaram, 2005;Nemoto et al., 2009). Thus, heterogeneous fracture transmissivity results from a combination of physical processes, such as shear deformation, orientation-dependent compression, and nonuniform preferential chemical alterations (Lander & Laubach, 2015;Lang et al., 2015).

10.1029/2018WR023189
Field observations of fractured rock permeability often reveal pronounced anisotropy, as a result of network geometry and state of stress (e.g., Barton & Quadros, 2015). Field studies have also shown that fracture-fault systems may develop higher permeability due to structural changes of the network, in the direction of fracture-fault intersection orthogonal to fault slip vectors (Sibson, 1996). Observations have led to empirical relationships that relate fracture hydraulic apertures to the physical process of shear displacement and the state of compression, that is, the acting normal stress (Barton et al., 1985). Hydromechanical numerical models using such empirical relationships have found that the permeability of fractured rock masses varies as a function of density and scale Min, Rutqvist, et al., 2004;Zhang & Sanderson, 2002;Zhang et al., 1996), as well as on stress-dependent mechanisms (Baghbanan & Jing, 2008). Remote stresses preferentially activate the deformation of specific fractures with a certain orientation to the far-field stress. These fractures can form high-conductivity pathways that span the entire domain and dominate the macroscopic permeability. Two-dimensional studies have found that the direction of maximum permeability of a fractured rock mass is aligned with the orientation of the maximum principal stress (Baghbanan & Jing, 2008;Jing et al., 2013).
The permeability of fractured rock masses is often modeled as a function of fracture network geometry (Bogdanov et al., 2003;Davy et al., 2006, De Dreuzy et al., 2001Ebigbo et al., 2016;Mourzenko et al., 2011;Olson et al., 2009;Saevik et al., 2013Saevik et al., , 2014. The transmissivity of individual fractures is often defined as a constant, is correlated with fracture length, or follows a direction-dependent distribution. Suggested universal relationships between flow and stiffness (Petrovitch et al., 2014;Pyrak-Nolte & Nolte, 2016) treat fractures as isotropic void spaces. Fracture transmissivities can also be defined as a mean or locally varying transmissivity field (De Dreuzy et al., 2012;Hamzehpour et al., 2009) or as a function of locally varying apertures resulting from the deformation of the rock mass in response to local or remote stresses (Salimzadeh et al., 2018). Therefore, the permeability of a geometrically isotropic fracture network may vary as a function of in situ stresses. This paper investigates the permeability of three-dimensional discrete fracture networks as a function of geometry and remote compressive stresses, with the aim of studying the effect of in situ stresses on the anisotropy of individual fracture transmissivity and its subsequent effect on the overall permeability of the fractured rock mass. The approach combines the simulation of processes at the mesoscale and at the small scale. Fractures are modeled as surfaces at the mesoscale, and contact stresses are resolved assuming linear elastic constitutive behavior of the rock. Then, fracture stiffness and anisotropic fracture transmissivities are computed locally for each individual fracture, as a function of contact stresses, using a small-scale model of the fracture. The method captures fracture roughness explicitly and assumes that the fracture is a self-affine composite surface that acts as an elastic body under confining pressure. The permeability tensor of the fracture-matrix system is then computed at the mesoscale, for several cases. Permeability of the system is shown to be highest in the direction aligned with the intermediate stresses, as a direct effect of the anisotropic fracture permeability that arises during contact.

Methodology
Discrete fracture models (DFMs) typically represent fractures as planar surfaces (Bogdanov et al., 2003;Koudina et al., 1998). The effects of small-scale features, such as roughness, are incorporated by means of parameters, such as transmissivity and friction coefficient. In the presented framework, a novel downscale-upscale approach is developed. Shear displacement and normal stress on fractures is computed at the DFM scale, in response to far-field stresses and the hypothesized opening-closure processes. Based on these values, individual fracture-scale models are created numerically to compute the transmissivity that (c) Trace of measured surface data of a Novaculite rock and (d) using radially averaged power spectra. The mismatch roughness wavenumber, q m , and the decreasing ratio between the spectra of the upper surface, C u , and composite surface, C c , illustrate the decreasing correlation between the upper and lower surface with decreasing roughness wavelength. (e) The discontinuous axis-averaged power spectra of the measured aperture fields, which are approximated in the model by (f ) smooth spectra for stochastic surface generation. results from the displacement and compression. Fracture transmissivities are fed back to the DFM, where global flow problems can then be solved. The full macro-scale permeability tensor of a rock mass can be computed for arbitrary fracture networks, with the friction coefficient as the only adjustable parameter. At this scale, the method accounts for flow through the matrix and fractures simultaneously and assumes that the matrix is permeable. As matrix permeability gets larger, the influence of the fractures on the overall permeability diminishes.
The method to compute the permeability of fractured rock masses has been validated by Lang et al. (2014). The validation of the rough surface generation against experimental data, as well as validations in the context of dissolution and precipitation, was presented by Lang et al. (2015). Validations in the context of fracture closure experiments are presented in Lang et al. (2016). The validation of the contact model has been presented by Nejati et al. (2016).

No Shear Displacement-Isotropic Aperture Fields
A fracturing process produces two isotropic, self-affine surfaces of near-Gaussian height distribution (Brown & Scholz, 1985;Poon et al., 1992;Schmittbuhl et al., 1995). The roughness power spectrum C(q) characterizes both the root-mean-squared roughness and the height correlation of such surfaces (Nayak, 1971). Spectra of this kind can be evaluated using a Discrete Fourier Transform and radial averaging of measured surfaces. A numerical algorithm based on a mathematical recipe in Persson et al. (2005, pp. 76-79) is implemented for this purpose. Spectra of this kind contain two parameters often used as characteristics for rough surfaces. The root-mean-square roughness, h rms , controls the normalization of the power spectrum (Persson et al., 2005). The Hurst exponent, H, reflects the surface height correlation length and corresponds to the slope of the power spectrum on a log-log plot. The Hurst exponent is related to the fractal dimension D f of the surface through H = 3−D f . It follows that H controls what may informally be described as the jaggedness of a surface. Further implementation details can be found in Lang et al. (2015Lang et al. ( , 2016. When brought into contact, without lateral displacement, that the two opposing surfaces should produce a perfect match and thus seal tightly. However, experiments in which tensile fractures are artificially induced, and the produced halves are reassembled to achieve maximum matedness, have always found significantly enhanced permeability (e.g., Gale, 1982;Moore et al., 1994;Morrow et al., 2001). This results from a mismatch between the two surfaces at small roughness wavelengths (Brown et al., 1986), mainly due to surface damage that seems to always occur during the fracturing process. It follows that while the opposing surfaces are correlated at the larger wavelengths, below some mismatch length scale, m , they are increasingly uncorrelated (Glover et al., 1998). This mismatch results in an aperture field that is responsible for the observed transmissivity of a tensile fracture.
The mismatch between surfaces of a tensile fracture can be illustrated using raw data (Yasuhara et al., 2006); see Figures 1a and 1b. After fracturing, these surfaces are cleaned, measured, and reassembled to maximum matedness. The aperture field of the fracture under negligible pressure, that is, when the surfaces touch at a single point, can be obtained by addition of the two profiles with opposite reference directions; see Figure 1a. The constructed surface is usually referred to as composite surface (Brown & Scholz, 1985) and is related to the aperture field under negligible elastic deformation, that is, for rigid surfaces through a = h max − h, where h max is the highest point of the composite surface and thus the point of singular contact. The main characteristic of a composite surface of tensile fractures is the lack of large wavelengths. This is a result of the strong correlation between the opposing surfaces. The aperture field shows a variable Hurst exponent: above Water Resources Research 10.1029/2018WR023189 a mismatch wavenumber q m , its scaling exponent is H = 0.8 and the spectrum of the composite surface coincides with those of the original surfaces ( Figure 1b). This indicates that above this frequency, the two surfaces are uncorrelated (Glover et al., 1998;Matsuki et al., 2006). Below this mismatch frequency, the surfaces are increasingly correlated, so that roughness at these wavelengths is not present in the aperture profile (see Figure 1e for measured and Figure 1f for modeled power spectra). In general, two random surfaces in contact, entirely uncorrelated, form a composite surface with identical Hurst exponent, for example, H = 0.8, with a reduced surface roughness amplitude h rms (Hyun et al., 2004).
The lack of large roughness features in the aperture profile of fractures without shear displacement has led to speculation that above some length scale, hydraulic and transport properties should remain invariant (Matsuki et al., 2006). Matsuki et al. (2006) devised a numerical approach complementary to the surface generation to create opposite surfaces that reflect the wavelength-dependent correlation between two tensile surfaces.

Shear Displacement-Anisotropic Aperture Fields
The mechanism responsible for mismatch between opposing surfaces in fractures without shear displacement is isotropic in nature. This is not the case in fractures with significant shear displacement, where a much larger mismatch results from relative displacement between the surfaces. The directional nature of this displacement results in an anisotropic void space (Figure 2), where flow channels tend to show more pronounced percolation perpendicular (⟂) to the direction of relative displacement than parallel (∥) to it. The ensuing anisotropy of fracture transmissivity has been measured as 1-1,000 (Gentier & Hopkins, 1997;Matsuki et al., 2006;Nemoto et al., 2009;Yeo et al., 1998) for laboratory-scale fractures. These experiments are set up to preserve the original surface roughness and to not allow for abrasion-related damage during the shear displacement. The assumption of surfaces remaining undamaged during the shearing process is also made throughout the present study. In natural rock, however, surface damage due to significant shearing may change both the measured prefactor in the h rms length scaling, and the scaling exponent H, in an anisotropic manner (Candela et al., 2012).
A statistical description of the shear-induced composite surface is difficult to achieve (e.g., Matsuki et al., 2006). Its anisotropic nature renders radially averaged power spectra ill suited for characterization. The established alternatives are linear power spectra of traces along the fracture surface and subsequent averaging of all trace Power Spectral Density of the same direction. This procedure results in separate C ⟂ (q) and C ∥ (q) spectra. The difficulty in stochastic generation of the aperture field results from a discontinuous spectrum in the shear-parallel direction; see Figure 1d. In general, the composite surface shows a roll-off similar to those in fractures without shear displacement, which reflects the fact that the opposing rock surfaces are correlated at wavelengths greater than the displacement magnitude.
An approximative method to generate composite surfaces similar to those seen in natural fractures is outlined by Brown (1995), using a single value for H and an anisotropy ratio between the x and y axes. Candela et al. (2009) devised a similar algorithm which directly incorporates H x and H y as input values but is restricted to fracture surfaces where no maximum roughness wavelength is defined. Both approaches have been implemented here, along with obtaining the composite surface directly from two mismatched surfaces and a given value for the shear displacement T ; see Figures 1c and 1d. This approach of generating two mismatched surfaces, displacing them with respect to each other, and bringing them into contact as a response to compression, combines methods presented by Matsuki et al. (2006) and Borri-Brunetto et al. (1999) and will be used herein to represent shear fractures and their contact .
Numerical solutions for the elastic contact between rough surfaces produce power law contact size distributions ( Figure 2c) similar to those of direct observations (Dieterich & Kilgore, 1996). Furthermore, aperture distributions under elastic compression ( Figure 2b) agree qualitatively well with those reported in experiments (Sharifzadeh et al., 2006).

Contact-Single Fracture
The frictionless contact between elastic bodies with rough surfaces is known to depend only on the aperture field before loading (e.g., Greenwood & Williamson, 1966). It follows that this contact problem can be reduced to the equivalent contact between a rigid realization of the composite surface ( Figure 1a) and an elastic, flat body; see Lang et al. (2016). As discussed, the rigid composite surface has a surface height profile that reproduces the original aperture field upon contact with the planar surface. The deformable body, correspondingly, has an elastic composite modulus E * , also referred to as reduced elastic modulus (Brown & Scholz, 1985), defined as E * = E∕2(1 − 2 ), where E and are the rock's native Young's modulus and Poisson's ratio, respectively. The numerical solution to the nonlinear contact problem implemented in this study is based on a fast Fourier transform (FFT) of the stress displacement integrals (Stanley & Kato, 1997). The contact pressure is transformed to the frequency domain by an FFT algorithm, multiplied by stiffness coefficients from the Westergaard (1939) solution for contact pressure due to a sinusoidal height profile, and the resulting displacement solution is brought back to the space domain by an inverse FFT. The inherent nonlinearity of the elastic contact problem that arises from the mutual dependence of contact area, contact pressure, and deformation is addressed within a conjugate gradient solver. This form of algorithm (e.g., Sainsot & Lubrecht, 2011) has been applied in numerous studies to solve the frictionless contact between rough surfaces (e.g., Almqvist et al., 2007;Jackson & Green, 2011;Lang et al., 2016;Sainsot & Lubrecht, 2011;Yastrebov et al., 2015).
The underlying grid is square, with  ×  rectangular cells of uniform size. In response to an external confining pressure, p, the obtained solution yields the contact stress, (x), and the normal displacement, over the fracture plane, x = (x, y). A mechanical aperture field, a m (x), is defined as the cellwise separation between the rigid and the elastic surface, with zero apertures at points of contact; see Figure 1b.

Flow-Single Fracture
Given the two-dimensional mechanical aperture field a m (x) as the local separation between the rigid and the elastic surface, a depth-averaged fluid pressure field, p f (x), is obtained by solving the Reynolds lubrication equation (e.g., Zimmerman et al., 1991): where f is fluid viscosity. As the numerical value of f is constant over the fracture, it has no influence on the steady state solution of the pressure field and can be set to unity. Any zero mechanical aperture of contacting elements is replaced by some small aperture value, w, representative of the permeability of the matrix rock.
Equation (1) is solved using a cell-centered finite-volume method. Two solutions, I and II, are obtained by imposing two macroscopic flow boundary conditions along the x and y axes, respectively. Specifically, the fluid pressure is prescribed to be constant along two opposing boundaries perpendicular to the main flow direction, and no-flow constraints are assigned to the flow-parallel boundaries.
The effective transmissivity of the fracture model is then obtained by Jaeger et al. (2007) and Tsang (1984): (1) calculating the flow rate vector field u(x) using Darcy's law as u = a 2 m ∕(12 f )∇p f ( Figure 2); (2) integration of the flow rate component along the macroscopic flow direction over the inlet cells to find the total volumetrix inflow Q, for example, over the left edge for flow parallel to the x axis (Figure 2), and (3) calculation of the fracture transmissivity, K, as where W is the width of the fracture inlet, or edge length of the quadratic composite surface model, and ∇p f is the macroscopic fluid pressure gradient. For fractures without shear displacement, K is a scalar property obtained as the harmonic average of K I and K II . K I and K II are the transmissivities of the fracture computed for two boundary conditions: flow along the x and y axes, respectively. For an isotropic system, these two values should be the same. However, the finite size of the system, and discretization errors, will inevitably induce small differences between these two permeabilities, and so the harmonic average of K I and K II is taken to represent the fracture transmissivity. For fractures with shear displacement, transmissivity is a diagonal tensor K with components K ⟂ and K ∥ , in the perpendicular and parallel direction to relative displacement, respectively. The corresponding parallel-plate hydraulic aperture, a h of the fracture in either direction can be obtained through the relationship (3)

Contact-Fracture Networks
The normal displacement, and, by extension, the normal stiffness between rough surfaces, is known to depend predominantly on the largest roughness wavelengths (Campañá et al., 2011). This allows for direct, (b) Subsequent triaxial compression leads to fracture closure under contact pressure p, shear stress , tangential displacement before contactĝ T , and tangential displacement in contact g T . Concurrently, individual models at the rough fracture scale (c) are constructed and subsequently (d) displaced and compressed based on network-scale average parameters per fracture. The confining pressure,p, is the harmonic average of the local pressure p over the surface. In (d) a m is the mechanical aperture due to elastic compression and is the discrete contact normal stress. S is the far-field stress tensor, and T is the shear displacement between fracture surfaces.
first-principle modeling of this problem by explicit representation of rough surfaces. As larger and larger wavenumbers are included in the surface representation, that is, q 1 is being increased, normal contact results converge, similar to results of transmissivity between rough surfaces (Neuville et al., 2011;Schmittbuhl et al., 2008). This is not the case for shear displacement, that is, frictional contact. Many attempts have been made to extend the normal contact approach to account for direct simulation of small-scale shear displacement (Bucher et al., 2002;Karpenko & Akay, 2001), but it is generally believed that there is no first-principle model for friction, which depends on multiple mechanisms that act across all roughness wavelengths (Persson, 2001(Persson, , 2005. It follows that friction is a phenomenological concept at the scale where the surfaces in contact are approximatively planar. It has been shown that at high levels of normal stress, the surface characteristics of rock do not seem to play an important role, and the friction coefficient converges to ≈ 0.6 (Byerlee, 1978).
The frictional contact at the DFM scale is solved to model the effects of far-field stress on the transmissivity of individual fractures and its effect on the permeability of a fractured rock mass. Fractures at this scale are modeled as split surfaces, which in numerical terms means that a duplicate set of nodes exist at fracture elements, one for each side of the fracture (Paluszny & Zimmerman, 2013). When modeling split fracture surfaces within the finite element framework, there is no explicit connectivity across the fracture, and contact between the two sides is accommodated by a supplementary problem based on an isoparametric contact discretization and a gap-based augmented Lagrangian algorithm. The formulation is described in detail in Nejati et al. (2016). The material behavior is assumed to be linear elastic, homogeneous, and isotropic. The finite element model makes use of an unstructured mesh based on quadratic tetrahedral elements to discretize the domain.
The surfaces of the model, including the internal fracture surfaces and external ones, are discretized using quadratic triangular elements. Conformal meshes are employed over the two surfaces of the fracture, leading to an element-to-element contact discretization based on triangular elements. The elements of the two fracture surfaces, which are often referred to as master and slave surfaces, are paired, which allows enforcement of contact constraints between each pair of elements. For quadratic models, the singular stress field near the crack front is modeled using quarter-point tetrahedra (Nejati et al., 2015a(Nejati et al., , 2015b. The frictional contact constraints are implemented in terms of an augmented Lagrangian-Uzawa method, details of which can be found in Nejati et al. (2016).

Hypothesized Opening-Closure Process
Contact between fracture surfaces on the DFM scale is computed by applying a set of sequential boundary conditions, resulting in a form of opening-closure history; see Figure 3. First, the initial opening of fractures is modeled by applying an initial fluid (normal) pressure p f over the fracture surfaces and solving for the ensuing deformation. An initial normal gap,ĝ N , between opposing surfaces, results ( Figure 3a). Second, the fluid pressure is set to zero, the solid deformation due to fluid pressure is discarded, and deformation and frictional contact is computed for an external triaxial compressive stress applied to the cube's boundaries. The far-field stress tensor is applied (Figure 3b), and S V > S H > S h is assumed. As a result of this compression, individual fracture surfaces may or may not go into contact. Moreover, during closure, crack surfaces experience a relative tangential displacement before and after coming into contact. The value of tangential displacement before contact is called the initial tangential gap,ĝ T , and any subsequent slip under compression is denoted by g T . The magnitude of the total local shear displacement between the surfaces is defined as T = | T | = |ĝ T + g T |. Furthermore, the normalized shear displacement ′ T is the dimensionless quantity for shear displacement, T . It follows that all dimensionless quantities are denoted by the prime superscript. If in contact, local values for contact pressure, p, and shear stress, , are then obtained (Figure 3b). Conceptually, this process is thought to mimic fracture initiation due to short-time, localized tensile stresses and subsequent closure due to the dominantly compressive state of stress that exists underground.

Permeability
To arrive at a full equivalent permeability tensor of a fractured rock mass, a series of single-phase, incompressible flow problems are solved over a porous medium domain using the finite element method. For each simulation, the resulting scalar pressure field is postprocessed to yield a pressure gradient vector field. Elementwise constant fluxes are computed using local element pressure gradients, and permeability is computed by applying Darcy's law. Contributions to the equivalent permeability tensor components are a function of element volume-weighted averaging of pressure gradients and fluxes. The resulting overconstrained problem is solved using a least squares approximation (Lang et al., 2014).

Sampling Domain
Averaged flux components are computed only over a restricted subvolume Ω of the flow region, away from the borders (Lang et al., 2014). This approach is referred to as subsampling, since the pressure diffusion problem is solved over the entire model, but the equivalent permeability is calculated by sampling the fluxes and pressure gradients in a subdomain only. Other authors refer to similar approaches as flow jacket or oversampling and find that this procedure reduces boundary effects, better reflects the connectivity of the sampling domain (e.g., Holden & Lia, 1992;Wu et al., 2002), and can be used to selectively evaluate the permeability of rock subdomains (Azizmohammadi & Matthäi, 2017). The sampling subdomain will be chosen such that only those elements are included that maintain a minimum distance to all model boundaries (Lang et al., 2014). Once the fluid flow problem is solved, the permeability of any number of subsampled regions can be evaluated. This allows investigation of homogenization scales, walking averages of permeability, and optimizing mesh refinement for the upscaling of permeability.

Two-Scale Coupling
Concurrently with the DFM solution to the frictional contact problem, a small-scale, rough surface model is constructed for each fracture in the network. This model consists of two fracture surfaces (Figure 3), one of which contains a mismatch factor as described in Lang et al. (2015). The rough surface depends only on the size of the fracture, which determines its roughness. Based on the harmonic average of shear displacement over the DFM fracture,̄T, the two surfaces are offset with respect to each other along the x axis, and the elastic, frictionless contact between them is computed under a confining pressure,p, which is also obtained as the harmonic average of the local pressure p over the surface; see Figure 3d. With the obtained mechanical  (Nemoto et al., 2009) at the same conditions. aperture field, the transmissivity of the fracture is defined as a fracture parameter K, which is a diagonal tensor of the form which is aligned with a local coordinate system of the fracture, where the x axis is parallel to the direction of the shear displacement vector T , the y axis is perpendicular to this vector on the fracture plane, and the local z axis is perpendicular to the fracture plane, that is, coaxial to the fracture plane normal vector. The finite element meshes for the mechanical deformation and fluid flow solution are identical, except for the fracture representation, which is split in the former and continuous in the latter. This is accomplished internally by merging split fracture walls.

Tensorial Transmissivity and Scalar Stiffness
To study anisotropic hydromechanical properties, three experiments have been modeled with the presented hydromechanical framework for single fractures (Bandis, 1980;Nemoto et al., 2009;Watanabe et al., 2008). All have been conducted on freshly induced fractures and with shear displacement under unstressed conditions so as to not damage the initial rock surface roughness. The length L of the numerical fracture models is chosen as the shorter edge length in the case of nonsquare experimental fractures. The Young's modulus E equals the reported values, that is, 66 GPa for slate (Bandis, 1980) and 29.4 GPa for Iidate granite (Hashida et al., 1993). Poisson's ratio has been chosen as = 0.15 for all models. It follows that two model geometries are created, one of L = 100 mm (Watanabe et al., 2008) and another of L = 50 mm (Bandis, 1980;Nemoto et al., 2009).
The first comparison of numerical with experimental results is made for the Watanabe et al. (2008) flow experiment. Three assemblies with increasing lateral offset between the surfaces are performed on an induced tensile fracture in Iinada granite. Resulting fracture permeabilities are evaluated perpendicular to the shear direction by means of permeameter experiments. For each offset T , a series of permeability values are then obtained for increasing levels of confining pressure; see Figure 4a. A follow-up experiment using the same rock type and setup is conducted later (Nemoto et al., 2009) to visually characterize areas of discrete contact that form due to increasing compression. This experiment solely focussed on the mechanical properties of the fracture.
The numerical model that was used to reproduce the observed permeabilities in Watanabe et al. (2008) is generated as follows. Two 100 × 100-mm surfaces with  ×  = 512 × 512 cells are created and displaced with respect to each other (Figures 3c and 3d) by distances equal to the shear offset used in the experiments, that is, T = 0, 1, and 3 mm, respectively. The composite surface, or aperture space, for each setup is then obtained by direct addition of the offset profiles. As the numerical contact algorithm works on a square model, any excess surface in the shear-perpendicular direction is discarded. From the shear displacements of 0, 1, and 3 mm, resulting surfaces are 100 × 100, 99 × 99, and 97 × 97 mm in size, respectively. Successively increasing confining pressures are then applied to obtain subsequent elastic contact solutions, and the hydraulic apertures of the resulting compressed aperture space are evaluated parallel, a h,∥ , and perpendicular, a h,⟂ , to the shear direction, along with the harmonic mean aperture of these two values, a h,× . Finally, the resulting values for a h,⟂ are compared to the experimentally obtained values. This procedure is repeated with varying values for the roughness prefactor in h rms = L H=0.8 (Candela et al., 2012) and mismatch wavelength m for the upper surface. A reasonable match (Figure 4a) is obtained with = 0.0023 and m = 55 μm. Note that the transmissivity of the fracture without displacement, that is, T = 0, depends mostly on m , and the transmissivity of the shear fractures depends mostly on the roughness prefactor . The permeability of a fracture, as reported in Watanabe et al. (2008), is related to its transmissivity by k = K∕a h .
The resulting mechanical aperture space for the T = 3-mm, 10-MPa case is illustrated in Figure 4b; a qualitative comparison of the contact area distribution measured experimentally (Nemoto et al., 2009) under identical conditions using pressure-sensitive film is shown in Figure 4c, for 15 statistically identical fracture realizations. Simulations and experimental data agree in that contacts of smaller size are consistently more frequent than larger area contacts. There is fair agreement between the two distributions, and this mechanical characteristic arises independently from the hydraulic matching procedure. Note that the nature of the distribution, a power law, and its approximate slope, is more significant than its magnitude, or integrated value, since the latter depends on the absolute value of contact area. It is expected that this absolute value is scale dependent and will depend on the visualization method used in the experiment.
The hydraulic properties of the T = 3-mm case have been further evaluated for pressures up to 50 MPa; see Figure 5a. Discrete pressure values have been chosen at log-spaced intervals, since this has been shown to be a useful functional relation between closure and pressure (Goodman, 1976) and thus facilitates higher resolution in regions of high curvature, that is, regions of strongly varying stiffness. The tangent (as opposed to secant) normal stiffness is computed numerically within the same framework, as the incremental change in normal stress divided by the change in the average normal separation of the fracture walls Lang et al. (2016).
The anisotropy, a h,⟂ ∕a h,∥ , is pronounced already at low pressures and increases with increasing compression p, up to a value of ≈ 195∕120 = 1.625. This implies a permeability ratio of (1.625) 2 and a ratio between the two transmissivity components of (1.625) 3 . The fracture specific normal stiffness of the T = 3-mm model shows conventional linear behavior (Figure 5b), which is the expected functional form of stiffness between randomly rough surfaces (Akarapu et al., 2011;Bandis et al., 1983;Berthoud & Baumberger, 1998;Hopkins et al., 1990;Swan, 1983). Note that the deviation of the last stiffness value from the linear trend is merely an artifact of the first-order finite difference approximation used for the end points of the closure curves. The linear scaling reflects that, with increasing p, existing contact areas grow and new contact areas form, so that distributions of stress and contact surface region remain constant if normalized to the area of discrete contact (Persson, 2006;Sun et al., 1985). It follows that the contact ratio R c varies linearly with p (Dieterich & Kilgore, 1996;Goodman, 1976) and that the distributions of contact stress and contact patch size scale linearly with p when normalized to the nominal area (Campañá et al., 2011). This linear scaling of contact ratio with pressure applies in the regime of low loads and elastic deformation (Carbone & Bottiglione, 2008;Greenwood & Williamson, 1966;Yastrebov et al., 2015).
The same surface parameters obtained for the permeability match against the Watanabe model have been used to independently reproduce the closure behavior of a fresh slate fracture with shear displacement of T ≈ 1 mm (Bandis, 1980); see Figure 5e. This experiment is chosen based on the prerequisites that fractures have been freshly induced and that the shear offset has been conducted under no load to preserve the original surface roughness. The agreement between the numerical and experimental closure results is excellent (Figure 5c).
The pressure dependence of the hydraulic properties of the slate fracture (Figure 5d) is qualitatively similar to that observed for the granite model Watanabe et al. (2008). The absolute value of the perpendicular-parallel

10.1029/2018WR023189
hydraulic aperture ratio, however, is smaller in magnitude, specifically a h,⟂ ∕a h,∥ ≈ 120∕80 = 1.5. This reduced anisotropy reflects the smaller relative shear displacement, ′ T = 1∕50 (Bandis, 1980) as opposed to 3∕100 (Watanabe et al., 2008). This highlights the fact that the ratio of absolute displacement to fracture size is an important indicator for hydraulic anisotropy. In other words, the ratio between shear offset and largest roughness wavelengths is a key parameter in the prediction of mechanical and hydraulic properties of a fracture. The agreement between numerically and experimentally obtained specific stiffness is similarly good (Figure 5e), a nontrivial conclusion given that the direct differentiation of experimental closure data is difficult due to the inherent problem of noise. This is of particular relevance in the context of mature versus new rock discontinuities, because, as opposed to new fractures, mature fractures and faults will generally be larger, due to intersection and coalescence, and have been subjected to larger and more complex shear deformation history.
The stiffness and transmissivity of two equivalent fractures of different sizes are also investigated. Two fracture models of L = 250 and 500 mm have been generated with constant cell size such that  = 512 and 1,024, respectively. An absolute offset of T ≈ 4 mm has been imposed on both models, and their specific stiffness is identical for all practical purposes; see Figure 5f. The equivalence of stiffness between these fractures, however, does not translate to the transmissivity of a fracture. This primarily reflects that stiffness/compliance represents the change in fracture void space with respect to pressure but not the absolute amount of the void space. Relationships between stiffness and transmissivity, therefore, are based on scaling laws that account for the size of the fracture in some way (Petrovitch et al., 2014;Pyrak-Nolte & Nolte, 2016). The established notion that stiffness alone cannot be used as a proxy for transmissivity is extended here by the observation that there exists no unique value of transmissivity anisotropy either (Figure 5f ). This is a nontrivial finding, in that the ratio a h,⟂ ∕a h,∥ does not depend on the void area in an absolute sense, and both stiffness and transmissivity anisotropy depend on size and offset only.
It can therefore be concluded that the transmissivity anisotropy of a single fracture increases both with relative shear displacement and increasing compression. The former controls the extent of flow channels that form due to the relative offset between the surfaces, and the latter reflects the fact that, as contact zones grow with increasing compression, they tend to percolate in the shear-perpendicular direction first, which follows in the sealing of shear-parallel flow at high compression. The normalized shear displacement ′ T is a physically more relevant parameter than the absolute value T .

Single and Conjugate Fractures Embedded in a Rock Matrix
Two models with rotating fracture geometries have been generated to illustrate the fundamental behavior of the coupled hydromechanical, multiscale modeling approach for DFMs; see Figure 6. A planar tilted hexagonal fracture surface has been chosen for these models; this allows for a precise discretization of the fractures as a finite element surface, avoiding the small discretization errors that would occur if circular fracture shapes were used.
Of the models depicted in Figures 6a and 6b, 13 instances have been generated with different angles of the symmetry axis with respect to the horizontal. The opening-closure process of the models has been simulated with an initial opening pressure p f = 0.9 S 3 and a subsequent triaxial far-field compression applied over the external model boundaries of The stress ratio is S 1 ∕S 3 = 3.1, which is theoretically the highest value that does not induce sliding in a critically oriented fracture for a friction coefficient of = 0.6; see Jaeger et al. (2007, p. 47).
A surface mesh box has been placed at the center of the model, aligned with the center of the symmetry axis in the x-z plane (denoted as Γ in Figures 6a and 6b). The purpose of this internal box is to force nodes to exist over its surfaces; it acts as a set of internal boundaries, on which the pressure boundary conditions of the flow problem are prescribed. This is done for the following two reasons: first, to have the fracture(s) percolate across the flow region, so that the obtained permeability reflects the fractured subsection of the model and second, to allow the sampling domain (denoted as Ω in Figures 6a and 6b) to rotate with the fracture geometry for all values of , so that the shape of the sampled fracture-matrix domain remains the same. This allows for direct comparison of the obtained equivalent values, with a constant fracture length and fracture-rock volume ratio over all configurations .
The edge length of the cubic finite element models (Figures 6a and 6b) is chosen to be 10 × 10 × 10 m, and the fracture length is L = 5 m. The domain is discretized by isoparametric tetrahedra. The underlying shape functions are quadratic for the elastic deformation-contact problem and linear for the flow problem. The maximum element size is defined as 10∕15 = 0.66 and adaptively decreased to 10∕40 = 0.25 within the flow region Γ; this flow domain has an edge length of 3×3×3 m. The elastic moduli of the rock matrix are E = 60 GPa and = 0.15. Upon solution of the mechanical opening-closure process, the hydraulic properties of the fracture(s) are evaluated, based on the FEM fracture average shear displacement,̄T, and average normal pressure,p. The opposing fracture surfaces are generated with H = 0.8, h rms = 0.01L H ,  = 8, 192, and a mismatch length scale of m = 500 μm. The upper wavelength limit due to sample size L equals the roll-off wavelength 0 . The absolute shear displacement is applied, and the composite surface is generated and downsampled to  = 512, using a bivariate spline approximation for the fracture-scale contact and flow analysis. The obtained transmissivity tensor is then assigned to the fracture in the finite element model by tensor rotation to align with the shear direction̄T, so that the parallel and perpendicular transmissivities align correctly with the shear direction. The matrix permeability is set to 1 mD. The subsequent pressure solution of the flow problem is obtained within Γ, and the permeability of the rock matrix ensemble is sampled over the subdomain Ω; see Figure 6.

Embedded Single Fracture
The hydromechanical results across the two scales are best illustrated for the single fracture case. At the network scale, the mechanical results of fracture confining pressure,p ′ , and shear displacement,̄T , are obtained from the contact solution, and the maximum rock-fracture permeability, k ′ max , results (see Figure 6d). At the single fracture scale, the elastic compression of the pore space, Δā ′ E , the contact ratio, R c , and the hydraulic aperture, a ′ h , result (see Figure 6e). The hydraulic apertures over the range of fracture orientations can be normalized with respect to their maximum value to obtain a ′ h , which represents the shear-perpendicular hydraulic aperture. The elastic compression of the fracture void space is the fractional reduction of the void space due to deformation, that is, Δa ′ E = ⟨a m ⟩∕⟨a⟩. These follow from the network-scale solution to the frictional contact problem, that is, the normal pressure, and the shear displacement; see Figure 6d. The contact parameters on the network scale for a fracture can be normalized using the fracture length for the amount of shear displacement,̄T ′ =̄T∕L, and the normal pressure can be normalized using the magnitude of S 1 , such thatp ′ =p∕S 1 .
A number of fracture orientations are examined by varying the fracture tilt angle . The effect of on fracture aperture and contact stress is investigated.
There is no shear displacement between the fracture surfaces for the two orientations of = 0 ∘ and 90 ∘ , since negligible shear stresses act on the surfaces. It follows that the composite surface, or fracture pore space, is identical for both configurations. The higher transmissivity at = 90 ∘ results from the lessened compression, as this configuration has the fracture normal aligned with the direction of the least principal stress, S 3 ; see Figures 6a-6c. At = 0 ∘ , the fracture is exposed to the largest possible compression, as it is oriented perpendicular to the direction of S 1 ; it follows that transmissivity reaches its lowest value. In both configurations, fracture transmissivity is a consequence of the m -induced surface mismatch and is thus isotropic. Between the two end-members of , shear displacement between the fracture surfaces reaches a maximum between 52 ∘ and 60 ∘ ; see Figure 6d. This range also represents the critical orientation of the fracture, that is, the fracture orientation most prone to friction-induced failure. However, no sliding between the surfaces is observed, as they stick together once in contact. This stick condition, g T = 0, is consistent with the stress ratio S 1 ∕S 3 = 3.1 for a friction coefficient of = 0.6. Thus, not only frictional failure but also the closing process before contact forms the largest offset between the surfaces in the range of ≈ 60 ∘ . The induced aperture reaches a maximum at ≈ 60 ∘ due to an optimal combination of shear and normal stress, which allows for maximum displacement during fracture closure and minimal compression in contact.
The uncompressed, rigid aperture of a fracture is proportional to its shear displacement. The hydraulic aperture of the fracture, however, reaches its maximum ( Figure 6e) at a larger angle, ≈ 60-67 ∘ . This results from the relatively large compression that the fracture experiences at its maximum shear displacement. The optimal ratio of the combined displacement and normal pressure occurs at a slightly steeper angle, which reflects the nonlinear drop inp ′ ; see Figure 6d and the reduced elastic compression of the void space, Δā ′ E Figure 7. The dependence of contact ratio, R c , on pressure, p, and shear displacement. (a) Numerical results for surfaces with  = 1, 024 and L = 50 mm. (b) Experimental results up to 60 MPa derived from pressure sensitive film images for a 50 × 50-mm granite fracture (Nemoto et al., 2009). in Figure 6e. The contact ratio between the surfaces, R c , has its maximum at the no-offset configuration at = 0 ∘ and decreases with increasing inclination of the fracture. This reflects the different combinations of offset and normal compression through all fracture orientations. A slight increase in R c at the = 90 ∘ , no-shear configuration is observed (Figure 6e), because of the large effect that shear displacement has on the contact area-load relationship (Figure 7). The permeability of the fractured rock region Ω (Figures 6a-6c), is proportional to the hydraulic aperture of the single fracture in a non-linear manner (Figure 6d). The contribution of the rock matrix acts to mitigate large gradients. Since a h,⟂ > a h,∥ , it follows that k ′ max (Figure 6d) is aligned with the y axis, except for the two no-offset configurations, in which the orientation either lies in the y direction or in the x/z direction, depending on whichever of the latter is numerically larger. The described opening-closure process qualitatively corroborates the established notion that critically oriented fractures are also the most likely to be hydraulically conductive (Barton et al., 1995).

Embedded Pair of Conjugate Fractures
The equivalent permeability of a rock with a pair of conjugate fractures in three dimensions is investigated. The angle between the fractures is fixed, and the pair is rotated with respect to the remote stresses to investigate the effect of their orientation on the permeability of the system. For the pair of conjugate fractures (Figure 6b), the direction of maximum permeability, k max , of the fractured rock region Ω is always along the y axis, since shear, and thus the displacement T , occurs only in the x-z plane. The more interesting question relates to the direction of the intermediate principal permeability, k med , which has to be located in the x-z plane. A prime example of such networks is bedding-plane-perpendicular fractures that form conjugate sets in limestone beds along the British Channel coast (e.g., Figure 7 in Belayneh et al., 2006). In such bedded sedimentary rocks, flow occurs primarily along the bedding plane, and the direction of largest permeability is of interest. For 13 instances of different inclinations of the symmetry axis with respect to the x-y plane, the hydromechanical results for the orientation of the intermediate permeability eigenvector in the bedding plane are shown in Figure 6c. The two end-member configurations, = 0 ∘ and 90 ∘ , have the expected orientation of k med along the symmetry line, since both fractures experience identical amounts of shear displacement and compression. In terms of magnitude, the bed-parallel permeability k med is higher at = 90 ∘ than at = 0 ∘ . This is due to larger shear and lessened compression when the symmetry axis is aligned with the maximum principal stress.
As increases from 0 ∘ to 90 ∘ , that is, the orientation of the principal stresses with respect to the fractures rotates. The direction of k med tends to align with the fracture with larger offset and lower compression. This is most pronounced around ≈ 30 ∘ , where one fracture is oriented perpendicular to the maximum principal stress, S 1 . In this configuration, the horizontally aligned fracture is offset-free and exposed to the largest possible compression, while the other fracture is near-critically oriented, with large offset and lessened confinement as a result. Beyond this inclination, the average flow direction moves gradually back toward the fracture pair symmetry axis. The largest permeability of the fracture-rock ensemble region Ω is obtained for a symmetry axis that is inclined 60 ∘ with respect to the S 1 -S 2 plane, very similar to the single fracture case. At this configuration, k med is almost 3 orders of magnitude larger than at the two no-offset orientations, as illustrated by the length of the corresponding eigenvectors in Figure 6c.

Fracture Networks
A series of 265 statistically identical fractured rock models (Figure 8a) have been generated to study the effect of stress on networks of fractures. The outer boundaries of the finite element model form a cubic domain with edge length of 40 × 40 × 40 m. The fracture region is restricted to the interior, such that all fracture centers are located at a distance equal or greater than 1.5 × R max from all outer boundaries, where R max is the upper limit of truncation in the fracture radii distribution. This distance to the outer model boundaries serves the purpose of restricting the effect of the fractures to the interior deformation of the model and to avoid interference with the applied stress boundary conditions. Eight fractures are generated in each model, following a power law size distribution with exponent −2.0 (Bonnet et al., 2001) and upper and lower fracture length truncation of 20.0 and 6.66 m, respectively. In total, 2,120 fractures are created, each with a different initially random but statistically identical small-scale roughness. Analogous to the single and conjugate fracture models, internal fluid pressure boundaries, Γ, are created in form of a centered cubic domain of 10 × 10 × 10 m to assess the fractured rock permeability in a percolating domain. Figure 8 illustrates the multiscale hydromechanical concept developed herein. After the solution of the deformation-contact problem at the network scale as a response to outer stress boundary conditions (Figure 8a), the averaged values of confining pressure and shear displacement are obtained for each fracture (Figure 8b). Based on these results, each fracture's composite surface is created and used to solve for rough surface elastic contact and ensuing transmissivity (Figures 8e  and 8f ). Equipped with the resulting transmissivity tensor and mechanical aperture, the network-scale flow is determined (Figure 8c). The permeability tensor,k, finally results (Figure 8d), here for the internal 8 × 8 × 8-m subdomain, Ω (Figure 8c) embedded in the flow domain, Γ (Figure 8b).
As illustrated using the single and conjugate fracture models, shear displacement and compression are functions of fracture orientation, which certainly holds true for mono disperse fracture populations. However, the amount of shear displacement is also a strong function of fracture size (Figure 8c), a fact that contributes to the motivation behind length-dependent fracture transmissivity relationships, which are commonly used for flow modeling in DFMs (Olson, 2003). The developed model qualitatively reproduces the observations that larger fractures, and those that are near-critically oriented, are stronger hydraulic conduits. The physical explanation for these differences is the larger pore space that results from large displacements (Figure 8b), low contact ratios, and low levels of compression for fractures close to the Coulomb failure regime (Figure 9c). In contrast, fractures away from this characteristic shear-normal are characterized by a small pore space and large contact ratios. An illustration of the dependence of aperture on fracture size and shear displacement is provided by comparing Figures 8b and 8c.
To isolate the effect of orientation with respect to the far-field stress from the effect of fracture size on transmissivity, length-normalized hydraulic apertures of all fractures in the generated networks are plotted as a function of the fracture poles; see Figure 9a. This normalized hydraulic aperture is defined as a ′ h = a h ∕L, where a h is the mean hydraulic aperture of the parallel and perpendicular component, and L is the fracture length. The isolated effect of fracture orientation on transmissivity ( Figure 9a) shows a preference for fractures oriented critically in the x-y plane, which is the plane of least and maximum principal stresses, S 3 and S 1 , respectively. This can be further illustrated by plotting a ′ h based on their projected stress state on a Mohr circle; see Figure 9e. The population of normalized hydraulic apertures spans 1 order of magnitude, which has even more pronounced consequences for the flow through these fractures, due to the cubic scaling of transmissivity.
Since fracture transmissivity is largest in the direction perpendicular to the shear displacement, in most cases, the preferred orientation of the overall k max of the fractured rock mass tends to be aligned with the y axis, that is, the direction of the intermediate principal stress, S 2 .
Not all fracture networks exhibit k max aligned with S 2 . In cases where single fractures dominate the orientation of k max , the permeability tends to be lower and generally reduced by a factor of around 1/3. This reduction is most pronounced when the direction of k max approaches the direction of S 3 ; see Figure 9b. An illustration of the way in which single fractures control the permeability of the rock mass is shown in Figure 8d. In this case, two connected fractures percolate in the horizontal direction (Figures 8b and 8c), resulting in k max being aligned with the x-y plane. Thus, a subset of fractures controls the overall flow in these cases. Despite the isotropic orientation of the fractures and the fact that these two fractures are at the same scale as the rock mass, the preferential direction of k max aligned with S 2 over all realizations is pronounced; see Figure 9b. The S 2 aligned orientation of k max is not only more frequent, but the permeabilities in this area are also larger in magnitude. This trend should be expected to be even stronger for larger fracture networks and for fracture sizes well below the size of Ω, which would move the model toward effective medium assumptions.
In this paper, individual fracture transmissivity is anisotropic and is the result of evaluating a small-scale micromechanical fracture contact model. The finding that k max is preferentially orientated along the intermediate stress S 2 is significant, in that it is observed in a three-dimensional model that accounts for the locally anisotropic transmissivity of individual fractures. Simulations in two dimensions, where fracture transmissivity is isotropic, have also found that the direction of k min is aligned with S 3 (Baghbanan & Jing, 2008; Figure 10. Distributions of (a) maximum permeability magnitude over 265 Discrete fracture model realizations, where k ′ max = k max ∕k matrix and (b) of the hydraulic aperture, a h , of all fractures therein. Note the implicit log scaling of the abscissa in (a). Min, Rutqvist, et al., 2004). In these two-dimensional studies, the direction of maximum permeability is found to be aligned with S 1 (Baghbanan & Jing, 2008;Jing et al., 2013). In the present three-dimensional study, where both S 2 and S 3 are taken into account, and where anisotropic fracture transmissivity is modeled, the direction of maximum permeability is found to be aligned with S 2 . These results are consistent with field observations reported by Sibson (1996) that over time, fracture-fault systems develop higher permeability in the direction of S 1 , orthogonal to fault slip vectors. This is due to the direct incorporation of the anisotropy effect a h,⟂ > a h,∥ . These effects should be expected to be even more pronounced for large offset fractures, beyond the framework of linear elasticity and small deformations, upon which the opening-closure concept developed herein is based. The population of obtained maximum permeability values, k max , has been fitted to a normal distribution; see Figure 10a. The obtained permeability values span over 3 orders of magnitude. This is typical of equivalent permeability results (Kirkby & Heinson, 2017), as opposed to effective permeability results (see Lang et al., 2014, for a discussion on the difference between effective and equivalent permeability). The latter are typically characterized by a much narrower standard deviation for large, statistically identical populations (Mourzenko et al., 2011). The population of obtained hydraulic apertures over all realizations has been fitted to a log-normal distribution; see Figure 10b.

Conclusions
A multiscale method has been presented to model the mechanical deformation of both matrix and fractures, to obtain anisotropic fracture transmissivities of individual fractures within three-dimensional fracture networks. Anisotropic transmissivities inform the computation of the permeability of these fracture networks and capture the effect of compressive stresses on the overall permeability of the fractured rock mass. The presented approach is based on the direct representation of the structural complexity that characterizes fractured rock across many scales, including power law size distributions and fractal geometries. The three-dimensional analysis developed herein includes the coupled direct simulation of individual fracture transmissivities as a result of network-scale shearing, taking into account the progressive deformation of the initially random small-scale fracture roughness. Results indicate that the shear displacement-induced anisotropic transmissivity of fractures leads to k max in the direction of S 2 . Thus, the direction of maximum permeability is most likely to be aligned with the direction of the intermediate principal stress. For isotropic fracture transmissivities, the presented findings support previous conclusions of k max ∥ S 1 (Baghbanan & Jing, 2008, Jing et al., 2013Min, Rutqvist, et al., 2004), due to a preferential mechanical activation of S 1 aligned fractures.