Mechanical Factors Controlling the Development of Orthogonal and Nested Fracture Network Geometries

Orthogonal fracture networks form an arrangement of open well-connected fractures which have perpendicular abutment angles and sometimes show topological relations by which fracture sets abut against each other, thus forming a nested network. Previous modelling studies have shown that orthogonal fractures may be caused by a local stress perturbation rather than a rotation in remote stresses. In this study, we expand on the implications of these local stress perturbations using a static finite element approach. The derived stress field is examined to assess the development of implemented microfractures. The results show that the continuous infill of fractures leads to a gradual decrease in the local tensile stresses and strain energies, and, therefore, results in the development of a saturated network, at which further fracture placement is inhibit. The geometry of this fully developed network is dependent on the remote effective stresses and partly on the material properties. Saturated networks range from: (1) a set of closely spaced parallel fractures; (2) a ladder-like geometry; and (3) an interconnected nested arrangement. Finally, we show that our modelling results at which we apply effective tension, are equivalent to having a uniformly distributed internal pore fluid pressure, when assuming static steady state conditions and no dynamic fluid behaviour.


Introduction
Orthogonal fracture networks are an arrangement of opening mode (mode I) fractures (Bai et al. 2002;Olson et al. 2007). These networks consist of two fracture sets perpendicular to each other. Several simplified mode I fracture network geometries inspired by pavement outcrops studies (Rives et al. 1992;Gross 1993;Ruf et al. 1998;Bai et al. 2002) are shown in Fig. 1. In mechanically layered rocks, opening mode fractures often form perpendicular to and terminate against mechanical interfaces (Bai and Pollard 2000;Bai et al. 2000;Schöpfer et al. 2011).
In the subsurface, orthogonal fractures, when remaining open under in-situ stresses, provide pathways for fluid flow (Bai and Pollard 2001) particularly in case of low matrix permeability. The opening of these Mode I fractures is controlled by the local stress field. Therefore, understanding effective permeabilities requires knowledge of the local stress field and it's effect on the local fracture behaviour (Bai et al. 2002;Gale et al. 2014;Hardebol et al. 2015).
In linear elastic fracture mechanics (LEFM), mode I fractures align parallel to ′ 1 and perpendicular to ′ 3 . Therefore, an orthogonal or nested fracture network geometry such as shown in Fig. 1b, c requires a rotation of the principal driving stresses. Numerical models developed by Bai et al. (2002) showed that when both in plane principal stresses ( ′ H , ′ h ) are in effective tension, the formation of orthogonal cross fractures can be explained by a local stress switch. These authors suggested that the formation of orthogonal fractures is driven by local deviations in the stress field and requires no regional stress rotation. They showed that this change in orientation occurs because of the reduction in effective stresses ( ′ h ) perpendicular to the fracture wall. The reduction of the effective tensile stress is called the stress shadowing effect (Fig. 2), and is comprehensively described in the literature (Bai et al. 2000a, b;Bai and Pollard 2000a, b;Schöpfer et al. 2011;Tang et al. 2008). Figure 2a shows how a mode I fracture opens and alters the local stress field, in effect to a remote tensile stress. Due to this opening, the local stress field ( �L x ) becomes perturbed, and normal stresses at the walls become zero. As two closely spaced layer-bounded fractures interact, the two regions of local stress perturbations interfere and the effective tensile stresses in the x-direction (driving stresses) reach 0.0 MPa (Fig. 2b, c). The spacing at which these effective tensile conditions reach 0.0 MPa is called the critical spacing (Bai and Pollard 2000a, b). Consequently, due to the absence of local tensile stresses (Fig. 2c), further fracture infill below this critical spacing becomes highly improbable. Hence, the two interacting fractures have reached a saturated state. Bai et al. (2002) showed that under bi-axial tension, applied tensile tangential stresses (parallel to the fracture wall) ( �L y ) were less effected by stress shadowing, resulting in a local stress switch ( ′L y > ′L x ). These authors argued that this effective stress switch could explain the occurrence of orthogonal fractures in between two mode I fractures. Other mechanical modelling studies (Olson et al. , 2009Olson 2007) used these bi-axial tensile stress conditions to systematically develop mode I orthogonal fracture networks, and argued that these network geometries could develop in one single deformation phase.
The aim of this study is to simulate the feasibility of the development of a saturated orthogonal fracture network under one constant remote stress regime, using a relation between LEFM and sub critical crack growth, which is used in several other modelling studies (Ko and Kemeny 2011;Olson 2007;Pollard and Segall 1987;Shen and Rinne 2007). We will focus on how an orthogonal fracture network shapes and is shaped by the regional applied effective tensile stresses, available strain energies and material conditions, and compare how the saturated geometries correlate to published field studies. Finally, we discuss whether mode Local stress perturbations as a result of LEFM and the concept of stress shadowing. It should be noted that tensional stresses are positive and that all layers have the same material properties. a Local stress perturbation due to the opening of single crack. b Stress shadowing effect due to two closely spaced opening fractures. c Two critically spaced opening fractures. Local stresses in between reach compression 1 3 I fracture networks can be present within the subsurface and we test and discuss the feasibility of replacing effective tension caused by fluid overpressures ( z > P f > H ⩾ h ) with tension assigned to model boundaries.

Model Setup
In our 3D mechanical simulations, we will mimic the effective tensile stress state by directly applying tensile stresses in both horizontal directions ( ′ H , ′ h ). The largest effective tensile stress is ′ h and corresponds to the smallest principle compressive stress plus the pore pressure ( The second horizontal stress is ′ H and corresponds to ( � H = 2 + P f ). Mimicking internal pressures with effective stresses or effective strains on the boundaries of models is a widely used approach (Olson et al. 2009;Philip et al. 2005;Rijken and Cooke 2001), and concerning numerical modelling, results in far less complicated mesh geometries. This implies that, the results presented depict static solutions of the applied boundary conditions and implemented fractures. Therefore, this approach cannot represent any transient effects, such as, pore pressure diffusion and the resulting changes in the local effective stress field. For our modelling study, we use the ABAQUS finite element software for solving the system of linear elastic equations. The free meshing algorithm in this modelling tool allows for the implementation of complicated fracture geometries, for instance, letting us change the aspect ratio ( L∕H ), spacing ( S∕H ) or orientation of the layer-bounded fractures. These features consist of fixed planes implemented as crack seams, with two walls as freely separable faces. The fractures thus have duplicate nodes, on both opposite fracture walls which allow for opening or slip, according to Linear Elastic Fracture Mechanics (LEFM) (Irwin 1957). The model domain has a linear elastic rock matrix which follows Hooke's law for 2D plane stress isotropic materials or 3D isotropic and these materials are defined by a constant elastic Young's modulus (E) and Poisson ratio (ν). The model domains are fixed in such a way that no rotation nor translation is allowed. This is done by fixing two bottom nodes in any of the free directions (x, y, z = fixed), and fixing the bottom surface in the z-direction to counter the overburden stress (Fig. 3a). The presented results follow the engineering sign convention so that tension is positive and compression is negative.
In our models tension ( ′ H , ′ h ) is assigned in both horizontal directions at the side walls, where ′ h and ′ H will represent the regional tensile stress in the x-and y-directions, respectively ( Fig. 3a, b). The vertical stress ( ′ z ) corresponds to the effective overburden stress (Fig. 3a, b). We assume that our model is placed at 1.0 km depth, assuming a typical layer cake type of geological basin. The tensile rock strength ( T 0 ) (Fig. 3b), is set equal to the maximum applied remote tensile stress ( � h =5.0 MPa) and will remain constant in all models. In the presented models we will change the Poisson ratio (ν), horizontal stress ratio ( parameters influence the local stress state and ultimately the saturated network geometry. Furthermore, a sensitivity analysis is done on the fracture spacing over height ratio ( S∕H ) and fracture aspect ratio ( L∕H ) to test how the fracture geometry alters the local stress field. The list of materials properties, assigned stresses and their respective ranges in magnitude are summarized in Fig. 3b.

Brittle Failure and Sources of Tensile Stresses in Geology
Mode I fracturing occurs when the tensile rock strength T 0 is overcome. This tensile rock strength can be used to construct a failure envelope in Mohr space. The failure envelope follows the Griffith and Mohr-Coulomb failure criterion for tensile and shear fracturing, respectively: where is the shear stress, n the normal stress, T 0 the tensional rock strength and the coefficient of friction.
In essence, all principal stresses in subsurface are compressive due to the overburden stress, Poisson stress and if present, superimposed tectonic stresses. Therefore, tensile stresses and the resulting mode I failure require additional internal forces or tectonic strains, such as local extension in a folding layer (Liu et al. 2016), exhumation and the resulting expansion (Fossen et al. 2007) or an alteration of softer and stiffer layers and the resulting local stress perturbations (Bourne 2003;Guo et al. 2017). Uniformly distributed pore fluid pressures are also regarded to be an important source of effective tension within rocks (Olson 1993;Philipp et al. 2013), where internal fluid pressures take part of the total stress that acts on the porous rock and leads to a less compressive, more tensile effective stress state. The effective ′ ij tensor describes this state of stress by: In case of a fully connected pore volume of the stratified rock column, pore pressures follow the hydrostatic depth curve. However, under a variety of geological scenarios, such as hydrocarbon generation in rocks with low matrix permeabilities or disequilibrium compaction (Ouraga et al. 2018;Tingay et al. 2009), pore fluids within rocks can reach significant overpressures (Emery 2016;Hantschel and Kauerauf 2009;Imber et al. 2014;Philipp et al. 2013). These overpressures form a possible explanation for effective tensile stresses in overburdened compressed rock layers (Gale et al. 2010;Germanovich and Astakhov 2004). In Fig. 4, this change in the effective stress state due to pore fluid pressures is depicted in Mohr space. The implemented fluid pressures can shift the Mohr circle potentially in to the tensile domain until tensile failure occurs ( � h = T 0 ).

From Local Effective Stresses to a Theoretical Sub Critical Fracture Propagation Rate
While the formation of mode I fractures driven by effective stresses that reach the tensile rock strength ( T 0 ) is understood ( Fig. 4), describing closely spaced mode I fractures is more difficult, since local tensile driving stresses in between a set of parallel mode I fractures get greatly reduced (Bai and Pollard 2000) (Fig. 2). This implies that the development of closely spaced mode I fractures and therefore orthogonal networks often occurs under sub critical stress conditions, since the local effective tensile stress never reach the theoretical rock strength. Therefore, in geoscience, rock failure under these sub critical stresses is often described using the empirical theory of Sub Critical Crack Growth (SCCG) (Atkinson 1984),which is based on the relation between the stress intensity ( K i ), fracture toughness ( K ic ) and measured fracture propagation rates. This fracture behaviour is often observed in laboratory experiments. For instance, double torsion experiments show fractures developing subcritically when the stress intensity factors fall below the critical value, i.e., K i < K ic (Nara et al. 2012; Shyam and Lara-Curzio 2006) Therefore, in geology, SCCG is considered to be an important mechanism for the development of complex fracture networks (Olson 2007).
Our study uses linear elastic fracture mechanics to incorporate SCCG in to our static models. This is done by calculating the stress intensity ( K i ) on the implemented micro cracks, using the modelled crack opening and the LEFM Pf

Mohr-Coulomb Curve
compression tension T0 Fig. 4 Regional stress state depicted in Mohr-Coulomb and the modelled effective regional stress state caused by fluid pressures 1 3 relations first described by Pollard and Segall (1987). These authors showed that the local stress intensity can be calculated using Eq. (3): where E is the Young's modulus of the matrix, the Poisson ratio, Y is a geometrical factor (set at 1 in this modelling study), a is half the crack length and C open is the maximum opening on a crack plane. The calculated stress intensity on each crack tip can then be related to a sub critical propagation rate using the stress corrosion cracking power-law (Atkinson 1984). This power-law relation is shown by Eq. (4) and is used in multiple lab an modelling studies describing SCCG (Atkinson 1984;Ko and Kemeny 2011;Nara et al. 2012;Olson et al. 2009;Shen and Rinne 2007): where v is theoretical sub critical growth rate (m/s), v max is the maximum crack propagation rate, K I(tip) is the stress intensity on the crack tip, K Ic is the mode I fracture toughness (material property) and n is the sub critical crack growth index. We assume that K Ic can be related to the tensile rock strength and an initial crack length ( K Ic = YT 0 √ a initial , where a initial = 0.0125 m), and the magnitudes of n and v max are set to 35 and 100 (m/s), respectively, which is similar to other modelling studies (Ko and Kemeny 2011;Olson et al. 2007Olson et al. , 2009Rijken 2005;Shen and Rinne 2007). The implemented SCCG approach states that low local tensile stresses result in a low fracture opening, which results in low stress intensities on the crack tip and therefore low calculated fracture propagation velocities. This implies that the reduction of effective stresses (Fig. 2) would result in a saturated network at which no implemented fractures will grow. In this study we will not focus on the impact of changing the sub critical crack growth index ( n ), which is known to have a significant impact on the fracture spacing in mode I fracture networks (Olson 2004(Olson , 2007. Another method of measuring the fracturing potential within an elastic material is using the strain energy density ( u ), which represents the elastic work stored in an elastically deformed body (Jaeger et al. 2007). This elastic energy is related to the applied principal stresses and elastic material parameters which is shown by Eq. (5): where u is the strain energy per unit volume, E is the Young's modulus, is the Poisson ratio and ′ L i are the local principal stress components. This equation states that minimising the available strain energy, will inhibit further fracture development and results in a saturated fractured material (Kemeny 1985).

A Step-by-Step Modelling Approach: Developing an Orthogonal Network Under Fixed Remote Stress Conditions
In the following, we take the computed 3D stress field, apply the relation between LEFM and SCCG and calculate local propagation velocities of implemented micro cracks, to depict the development a saturated fracture network under fixed remote stress conditions, assuming steady state conditions and no transient pressure effects. For this modelling scenario, a remote tensile stress ratio of ( and a Poisson ratio (ν) of 0.25 are applied to the boundaries and matrix of the model, respectively. The results depict a z-slice through the model, essentially cutting the model in half.
The initial fracture network (Fig. 5) depicts two initially widely spaced parallel fractures ( S∕H = 4.0, spacing = 0.4 m) placed perpendicular to the highest effective tensile stress. In between these fractures, 17 microfractures are placed perpendicular to one of the principal stress directions. These microfractures have an initial length of 0.025 m (2*element size) and have an initial spacing of x = 0.1 m for parallel fractures, and y = 0.15 m for orthogonal fractures, respectively. The fractures are illustrated by the aligned circles, which size is proportional to the local fracture aperture. This geometry (Fig. 5) is used as an initial fracture network for all modelling scenarios, which aim is not to represent a naturally imperfect rock volume, but to model and depict the mechanics controlling the development of an orthogonal fracture network. In our numerical analysis, the two outer fractures are not allowed to propagate throughout each timestep and we will only focus on the development of fractures within the area of interest (Fig. 6).The initial network geometry and the highest effective tensile stresses are shown by Fig. 6a. This figure depicts how the implemented mode I cracks open and perturb the local stress field following LEFM. Using Eqs. (3) and (4), the theoretical propagation velocity (m/s) of each micro crack is calculated. As shown by the local stress field, local tensile stresses, and therefore, local fracture opening is highest at the centre of the model (x = 0.0, − 0.2 < y < 0.2). This implies that the parallel fractures in these locations are most likely to grow. The resulting network derived from the first timestep is depicted in Fig. 6b and shows that the fractures in the centre of the model have indeed developed into a fully grown parallel fracture. Furthermore, some other micro cracks showed a few cm of growth (Fig. 6b) ) are becoming shadowed following LEFM. However, the local principal tensile stresses in the y-direction remain present (Fig. 6b), resulting in a local stress flip ( ). Despite the fact that a local stress flip occurs, the network geometry calculated during this stage is still a set off closely spaced parallel fractures (Fig. 6c), which can be explained by the fracture growth implemented during the previous timestep.
At the third timestep, the newly implemented fractures further perturb the local stress field (Fig. 6c). This figure shows how the effective stresses in the x-direction are now completely shadowed ( �L x T 0 ≈ 0 ), inhibiting further parallel fracturing. However, as was observed during the previous timestep, tensile stresses in the y-direction remain present (Fig. 6c), causing the orthogonal microfractures to remain open. Given enough time (Myr) and prolonging tensile stresses, these fractures will develop into orthogonal fractures and this is also derived by our model (Fig. 6d).
The results depicted by the final stage (Fig. 6d) show a saturated (fully developed) fracture network, where all initial tensile stresses are shadowed as a result of the emplaced fractures, essentially implying that: ≈ 0. This local stress field, therefore, indicates that no new mode I fracture can grow in between the current active fractures (as long as these fractures remain open and perturb the stress field). Hence, this figure represents the final network geometry which is most probable to form, regarding this specific initial network geometry, remote stress state, implemented Poisson ratio and sub critical crack growth parameters.
The development of a saturated orthogonal network (Fig. 6) also results in a minimization of the local available strain energies (Eq. 5). The step-by-step reduction in available energy as a function of continued mode I fracture infill is shown in Fig. 7. At the final development stage (Fig. 7d), the strain energy density created by the applied tensile stresses is reduced to 0, inhibiting further mode I fracture infill. Therefore, Fig. 7d depicts a mode I fracture network geometry at which the available strain energy density is minimized, and, therefore, has the highest probability of forming under these initial conditions.

Sensitivity Analysis
The development of a saturated orthogonal network is a function of specific boundary conditions (Figs. 6, 7). Therefore, it is important to apprehend how changing these boundary conditions, alters the local stress perturbations and fracture behaviour. In the following, a sensitivity analysis is performed which tests how altering the applied boundary conditions (interacting fracture spacing: S∕H , applied tensile stresses:

The Effect of the Applied Tensile Stress Ratio on Local Stresses
The ratio of the applied horizontal stresses (  Figure 8a shows the resulting stress spacing curve which depicts the local measured stress as a function of the ratio between the applied stresses and the relative fracture spacing. As was shown by Bai and Pollard (2000), tensile stresses in the x-direction [ ′L x ], occur until a fracture spacing of S/H = 1.0, lower spacing ratios show a switch to compression [ �L x < 0.0 MPa] . The applied tensile stresses in the y-direction are less affected by this stress shadowing effect, hence a local stress switch occurs ( ′L x < ′L y ). As expected, the relative fracture spacing and at what stress magnitude this stress switch occurs is highly dependent on the applied tensile stress in the y-direction ( ′ H ) (Fig. 8a).

The Effect of Poisson Ratio on Local Stresses
Linear elastic fracture mechanics indicates that tensile fracture opening results in zero normal stresses [ �L x = 0 ] and positive local strains [ L x > 0 ] (opening of the fracture) at the fracture walls. Therefore, following from the tensor behaviour of the compliance matrix, fracture opening will result in locally generated compressional tangential stresses �L y < 0] . The magnitude of this local stress response is dependent on the Poisson ratio of the matrix, with a higher Poisson ratio resulting higher magnitude of the response. Figure 8b shows this local stress results as a function of fracture spacing and the modelled Poisson ratio. A low ratio of 0.1 results in almost no local tensile stress loss in the y-direction, whereas a high Poisson ratio results in the opposite effect. Figure 8c illustrates the differences in local stress magnitude [ �L y , �L x ] with respect to the aspect ratio of the fractures. This figure shows that the aspect ratio of the layer-bounded fractures only has a minor effect on the local stress state.

3
This minor effect can be explained using the total fracture opening, which is higher for fractures with a higher aspect ratio. This increased opening results increased strain away from the fracture wall, which results in an increased compressional tangential stress response and, therefore, reduced effective tensional stresses. An increased opening also results in an increased area of stress shadowing (Pollard and Segall 1987) hence a lower [ �L x ] for fractures with a larger aspect ratio (Fig. 8c).

Modelling Results: Saturated Network Geometries as a Function of the Remote Stress State and Poisson Ratio
In the following, we use the same approach as described in Sect. 3.2, in order to investigate and depict how different boundary conditions lead to different saturated network geometries. The different modelling scenarios are shown in Table 1.
The first three scenarios test the effect of Poisson ratio on the development of an orthogonal fracture network ( Fig. 9a-c). Scenario 1 tests the effects of a high Poisson ratio on the saturated network geometry, the applied stress ratio is the same as shown in Sect. 3.2. As shown in Fig. 8b, a high Poisson's ratio, results in a high tangential stress response. Therefore, following our modelling approach, orthogonal fractures are unlikely to develop, and the modelled saturated network geometry resembles a set of critically spaced parallel fractures (Fig. 9a). In modelling scenario 2 a more common Poisson's ratio of 0.25 is implemented and this lower ratio results in more favourable conditions regarding the presence of local tangential tensile stresses (Fig. 8b) and, therefore, the development of the orthogonal microfractures. Hence, the saturated network geometry resembles a critically spaced set of systematic fractures (parallel) with orthogonal fractures in between (Fig. 9c). Scenario 3 shows the end-member results of a low Poisson's ratio. Following from Hooke's law materials with a low Poisson's ratio show a low tangential stress response (Fig. 8). Hence, orthogonal fractures are more likely to develop concerning our modelling approach. The final saturated network geometry, therefore, 1 3 shows a developed ladder-like fracture network geometry (Fig. 9c).
Scenarios 4, 5 and 6 test the effect of the applied remote stress ratio on the development of a mode I fracture network ( Fig. 9d-f). As expected due to the low remote stress ratio ( ′ H / ′ h = 0.4), orthogonal failure becomes highly improbable. Therefore, the saturated network consists of a set of critically spaced parallel fractures (Fig. 9d). Scenario 6 ( Fig. 9f) has isotropic horizontal stress conditions ( ′ H / ′ h = 1.0) and, therefore, resembles a favourable concerning the development of an orthogonal network. Due to the magnitude of ′ H , the first orthogonal fractures already develop after the initial stage, resulting in long orthogonal fractures in between the two initial fractures. These orthogonal fractures create a new stress shadow in the y-direction. However, local tensile stresses in x-direction are not completely shadowed; therefore, the development of the smaller implemented parallel fractures is not halted. In conclusion, the fully saturated network geometry resembles a nested orthogonal fracture network (Fig. 9f).

Modelling Effective Tensile Stresses: Can Boundary Tensile Stresses Replace Internal Fluid Pressure?
Applying effective tension on the model boundaries is a widely used method to simulate mode I fractures in the subsurface (Bai et al. 2002;Germanovich and Astakhov 2004;Olson 2007;Olson et al. 2009). However, effective tensile stresses in the subsurface, and therefore, mode I fracture development is often supported by pore fluid pressures (Olson 1993). Therefore, to test whether it is representative to replace this effective geological stress state with applied tensile stresses, we have created two 2D models. Both approaches are shown in Fig. 10.    Figure 10b shows the alternative, commonly applied model design with effective tensile stresses assigned to the model boundaries. Both models tested assume steady state conditions and represent the crack tip as a singularity point. This required a structured meshing approach, which is different from the free meshing approach used for the 3D models. The results (Fig. 10c-f) indicate that both approaches lead to an identical local stress distribution, which is expected from the superposition principle (Jaeger et al. 2007). This confirms that the application of effective tension on the boundaries is a valid representation of modelling the effects of homogeneous distributed fluid overpressures within the matrix, assuming steady state static conditions. However, our pore pressure models, and therefore the results presented in Figs. 6,7,8,and 9, cannot account for any dynamic pressure behaviour, such as pore volume change, heterogeneous pressure fields or any fluid diffusion due to the opening of the cracks. All of these factors play an important in role in the formation and development of natural hydraulic fractures (Mandl 2005;Ouraga et al. 2018;Philipp et al. 2013). Recent modelling work done by Ouraga et al. (2018) has shown that this dynamic pore pressure behaviour together with hydro-mechanical coupling results in transient fracture behaviour, at which fractures open and close as a result of local dynamic pressure changes. They showed that this dynamic fracture behaviour could have a significant impact on the resulting fracture network. Furthermore, their models showed more dynamic features such as small spacing/height ratios and non-elliptical fracture opening, which are also observed within natural fracture networks (Gudmundsson et al. 2010;Phillip 2008). However, these authors also showed that, given enough time, this transient behaviour ultimately resulted in similar network geometries and local effective stresses as depicted in this study and other mechanical modelling work (Bai et al. , b, 2002Olson et al. 2009;Ouraga et al. 2018). Nonetheless, it is important to understand, that when abruptly halted, transient fracture behaviour can result in different network geometries that the expected geometries derived from mechanical modelling.

Linking Network Geometries to Remote Stress, Local Perturbations and Material Conditions
Under steady state conditions and following the stress shadowing effect, the development of networks with closely spaced orthogonal fractures must occur under subcritical conditions (Bai et al. 2002;Olson 2007). Stresses in between two or more interacting fractures are greatly reduced and are, therefore, subcritical, with local conditions never reaching the theoretical tensile rock strength. This implies that without subcritical growth, no new fractures would grow in between two interacting opening mode fractures. However, when fracturing under sub critical stresses is acknowledged, the coeval development of orthogonal fractures in regions subjected to a local stress switch can occur, which can lead to the development of an intertwined nested network (Olson 2007). Our models used a relation between local stresses, fracture opening and propagation velocities to assess fracturing under subcritical conditions. The resulting models show how the Mode I saturated network geometry (parallel, ladder-like or nested) is highly dependent on the applied stress conditions and material parameters (Fig. 9). Our results showed that the successive development of an orthogonal fracture network requires sufficient remote tension in both horizontal directions and a favourable Poisson ratio ( ′ H / ′ h ≥ ± 0.7, ν ≤ 0.25). If these remote stress conditions are not satisfied, the saturated network geometry will probably be resembled out of a set of critically spaced parallel fractures. Our models also show that a small contrast between the applied horizontal stresses ( ′ H / ′ h ≥ 0.8) will most likely result in the development of a nested network; a network pattern similar to hierarchies in orthogonal networks observed in outcrops (Engelder and Peacock 2001;Mandl 2005;Olson 2008). These model findings help distinguishing parameter space domains (remote stress state and Poisson ratio) (Fig. 11) at which the development of one of the saturated network geometries is most probable.

Implications for Field Geology and Subsurface Reservoirs
The modelling results of this study show the conditions under which an orthogonal fracture network is most probable to form while experiencing one constant remote stress regime. The question is in how far these learnings can provide basic rules for predicting the presence of an orthogonal fracture network in the subsurface or on outcrops. First of all, we should acknowledge that in nature, the development of mode I fracture networks is driven by extensional strain rates instead of applied tensile stresses. Studies modelling network development as a result of these strain show complex orthogonal networks (Olson 2007;Olson et al. 2009), whereas our models depict relatively simple configurations. However, our networks do show how continued fracturing results in saturated network geometries at which further fracture infill is prohibited, and how this geometry is dependent on the applied boundary conditions. As was explained in Sect. 2.2, tensile stresses can be contributed to a variety of sources. Therefore, one could argue and assign no internal pore pressure. c Local horizontal stress distribution for our pore pressure model d Local vertical stress distribution for our pore pressure model. e, f Local stress distributions regarding our effective tensile stress models that orthogonal fracture networks, observed on outcrops, can be created by uplift and erosion, and the resulting extensional strains in both horizontal directions (Fossen et al. 2007). This implies that orthogonal fracture networks are erosional features and are, therefore, unlikely to be present in the subsurface. Furthermore, exhumation of rocks occurs at geological time scales, resulting in sufficient time for fracture development under sub critical conditions. Although the above stated argument is true, several recent field studies suggest that orthogonal fracture networks can also form as a result of fluid overpressures. For instance, recent studies on the low permeable layers of the Whitby Mudstone Formation and Cleveland Ironstone Formation (Cleveland basin, Yorkshire, UK) (Emery 2016;Imber et al. 2014), argued that the observed orthogonal and nested fracture networks were created by pore fluid overpressures. Furthermore, micro seismicity and core data from the low permeability reservoirs found evidence of open parallel and orthogonal fracture networks being present within the subsurface (Fisher et al. 2005;Gale et al. 2010;Warpinski 2009). The observed fracture networks were believed to be created by fluid overpressures prolonging for a relatively long time, resulting in effective tensile stress and the development of closely spaced orthogonal fracture networks (Fisher et al. 2005;Warpinski et al. 2014;Warpinski 2009). These studies also argued that the stimulation, effective permeability and resulting production was significantly enhanced by the presence of the natural fractures. Therefore, the mechanics and local driving stresses controlling the development of parallel, orthogonal and/or nested fracture networks, can for instance, be importing regarding well optimizations and hydraulic stimulations plans concerning tight reservoirs. However, our modelling results can not account for any pressure diffusion, and, therefore, should not be used for fracture network characterization in permeable reservoirs studies.

Conclusions
This study shows how orthogonal and nested network arrangements can originate from stress shadowing and the resulting switching of the maximum tensile stress, caused by closely spaced Mode I fractures. First, the systematic development of an orthogonal and nested network requires that both principle horizontal stresses can act as a tensile driving stress and, therefore, are in sufficient tension. Second, due to the opening of systematic fractures and the resulting reduction in local effective stresses, the simultaneous development of an orthogonal network occurs under sub critical stress conditions. Furthermore, our modelling depicts the underlying process on how the coeval development of mode I fractures (parallel and/or orthogonal) leads to a saturated fully developed network which can have a nested, orthogonal or parallel arrangement. As the network develops, local stresses and consequently local failure probabilities are greatly reduced by stress shadowing to the extent that eventually no new fracture can develop and a saturation of the network is reached.
Finally, this study confirms the mechanical feasibility of developing parallel, orthogonal or nested fracture network topologies under bi-axial effective tensile stresses caused by, for example, pore fluid overpressures. The models illustrate that the geometry of the developed saturated network is highly dependent on the applied regional stress state and Poisson Ratio of the material, when assuming a fixed initial fracture network. A small contrast in magnitude of the horizontal effective tensile stresses ( . The oblique angle of these different domains is based on our models which tested the effect of the Poisson ratio (Figs. 8,9) Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.