Conditions for fracture arrest in layered rock sequences

Fracture arrest in layered rock sequences is important in many geodynamic processes, such as dyke-fed volcanic eruptions, earthquake ruptures, landslides, and the evolution of plate boundaries. Yet it remains poorly understood. For example, we do not fully understand the conditions for dyke arrest (preventing potential eruptions) or hydraulic-fracture arrest in gas shales (preventing potential aquifer pollution). Here we present new numerical results on the conditions for arrest of fluid-driven (mode-I) vertical fractures in layered rock sequences when the tips of the fractures approach the interface between two layers of contrasting mechanical properties. In particular, we explore the stress-field effects of variations in layer stiffness, proximity of fracture tip to layer interface, and layer thickness. When the layer hosting the fracture tip is stiffer, fracture arrest normally occurs at the interface with the more compliant layer. By contrast, when the layer above the interface is stiffer, fracture arrest may occur within the host layer well below the interface. These conclusions are supported by field observations of arrested fluid-driven joints and dykes and, therefore, provide a better understanding of the mechanical conditions for


a b s t r a c t
Fracture arrest in layered rock sequences is important in many geodynamic processes, such as dyke-fed volcanic eruptions, earthquake ruptures, landslides, and the evolution of plate boundaries. Yet it remains poorly understood. For example, we do not fully understand the conditions for dyke arrest (preventing potential eruptions) or hydraulic-fracture arrest in gas shales (preventing potential aquifer pollution). Here we present new numerical results on the conditions for arrest of fluid-driven (mode-I) vertical fractures in layered rock sequences when the tips of the fractures approach the interface between two layers of contrasting mechanical properties. In particular, we explore the stress-field effects of variations in layer stiffness, proximity of fracture tip to layer interface, and layer thickness. When the layer hosting the fracture tip is stiffer, fracture arrest normally occurs at the interface with the more compliant layer. By contrast, when the layer above the interface is stiffer, fracture arrest may occur within the host layer well below the interface. These conclusions are supported by field observations of arrested fluid-driven joints and dykes and, therefore, provide a better understanding of the mechanical conditions for dyke-fed eruptions.

Introduction
Understanding and forecasting the conditions leading to propagation or arrest of fractures in crustal segments composed of contrasting layers is a fundamental unsolved problem in Earth Sciences. The solution of this problem is of great importance because fractures largely control key processes in fields such as volcanology, seismology, engineering geology, hydrogeology, and petroleum geology. Most rock fracture propagation models assume that the crustal segment hosting the fracture is either isotropic and homogeneous (i.e. non-layered), or, at most, composed of only a few layers or units ( Segall, 2010 ). For most volcanic eruptions to occur, however, a fluid-driven fracture, a dyke, must be able to propagate through numerous crustal layers and interfaces from its source magma chamber to the surface ( Gudmundsson, 2016( Gudmundsson, , 2020. Thus, the mechanical conditions that allow a propagating dyke to successfully penetrate all the layers ahead of it, rather than become arrested, provide one of the main controls on whether an eruption occurs or not. Also, man-made hydraulic (fluid-driven) fractures are commonly used to increase the permeability of fluid-filled reservoirs of various types. Such fractures propagate in a manner that is mechanically analogous to that of dykes, particularly the vertically oriented hydraulic fractures that are commonly generated to increase permeability in un- Here the thickness (palaeo-aperture) of the dyke is 0.25 m at the bottom of the exposure, and it thins towards the tip, where the dyke was arrested before reaching the interface with the stiff lava. (B) Series of joints in a interbedded limestone (Lst) (stiff) and shale (compliant) sequence at Nash Point, South Wales, UK. The thicker shale bed is approximately 0.3 m. Within the centre of that shale bed one of the joints is arrested before reaching the stiffer limestone above, whereas other joints have propagated through the whole sequence or become arrested at the interface between the two layers. (C) and (D) are schematics of (A) and (B), respectively, where black lines depict interfaces between different rock units and red lines depict the outline of the dyke in (C) and the joints in (D). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 2017 ; Montgomery et al., 2010 ;Yao, 2012 ;Virgo et al., 2014 ;Fisher and Warpinski, 2012 ;Kavanagh et al., 2015 ;Tang et al., 2018 ;Gudmundsson et al., 2010 ;Barnett and Gudmundsson, 2014 ), both using analytical ( Kavanagh et al., 2015 ;Kavanagh et al., 2017 ;Kavanagh et al., 2013 ) and numerical modelling techniques ( Virgo et al., 2014 ;Ghani et al., 2015 ;Zhao et al., 2017 ), the conditions controlling fracture arrest have not received nearly as much attention. Furthermore, many studies on fracture propagation do not fully represent the heterogeneous (i.e. layered) nature of the crust, even though many crustal segments are known to be heterogeneous.
The studies that do consider the influence of crustal heterogeneity include data taken from field observations, laboratory experiments, and analytical and numerical models ( Kavanagh et al., 2015 ;Barnett and Gudmundsson, 2014 ;Teufel and Clark, 1984 ;Maccaferri et al., 2011 ;Bonafede and Rivalta, 1999 ;Maccaferri et al., 2010 ;Taisne et al., 2011 ;Taisne and Jaupart, 2009 ;Kavanagh et al., 2013 ;Gudmundsson et al., 2002 ;Geshi et al., 2012 ;Brenner and Gudmundsson, 2004 ;Larsen and Gudmundsson, 2010 ;McGinnis et al., 2017 ;Smart et al., 2014 ;Douma et al., 2019 ) and they generally suggest that such heterogeneity does affect fracture propagation. In volcanology, some authors attribute dyke arrest to either a negative buoyancy contrast between the magma and the host rock, or an insufficient supply of magma from the original reservoir ( Maccaferri et al., 2011 ;Taisne et al., 2011 ). In a more general sense (but which also includes volcanology), most studies that consider crustal heterogeneity find that fractures become arrested at contacts, or interfaces, between layers of contrasting mechanical and elastic properties, regardless of buoyancy effects ( Gudmundsson, 2011 ). Generally, there are three contact mechanisms that influence fracture arrest which are: Cook-Gordon debonding (delamination), stress barrier and elastic mismatch. A complete explanation of each mechanism is provided in the Supplementary Information and can also be found in Gudmundsson ( 2011Gudmundsson ( , 2020. However, there have been numerous field observations which indicate that mode-I fractures can also become arrested before reaching the interface or contact between different layers ( Fig. 1 ). The above three mechanisms may then not be sufficient to fully explain these observations in which case new models are needed. The focus of this study is, hence, to understand how fractures can become arrested prior to meeting an interface between heterogeneous rock layers or units.
The example illustrated in Fig. 1 A is a basaltic dyke on the island of Tenerife, Canary Islands, whose palaeo-aperture (thickness) is approximately 0.25 m at its base and thins towards the tip. The stiffness (Young's modulus) of the tuff layer hosting the dyke is much lower than that of the lava flow above the tuff. In this section, the dyke became The model is 300 × 300 m in vertical cross-section -made so large to avoid any edge effects on the results. The 10 m long (dip dimension) fracture is fully contained within the middle layer (host layer) which has the same properties as the grey top and bottom layers. The properties of the two green layers vary throughout this study. (B) Zoomed-in image with annotations demonstrating which layer hosts the fracture and which is ahead (green). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) arrested prior to reaching the interface between the tuff and the lava flow.
The fractures illustrated in Fig. 1 B occur in an interbedded limestone and shale sequence at Nash Point, South Wales. The fracture pattern is ubiquitous through the entire sequence and along the several kilometres of coastal exposure, and consists primarily of joints (that may or may not be fluid driven) and mineral veins which were both formed as predominantly mode-I fractures. In this example, the limestone layers are much stiffer (with a higher Young's modulus) than the shale layers. In Fig. 1 B three fractures propagate across the sequence. Two of the fractures cross quasi-vertically through the entire sequence, while one becomes arrested within a shale bed, and does not reach the interface with the limestone layer above.
The effect of layering on tensile stress distributions at the fracture tip prior to meeting an interface between layers has received comparatively little attention. Müller (1986) concluded that the stress intensity factor at a fracture tip decreases when a fracture propagates towards a stiffer layer, but increases when the fracture propagates towards a more compliant layer. However, the effect of layer thickness was not considered, and some other important model parameters were omitted (such as fracture aspect ratio, internal fluid pressure, and the model dimension) which makes it difficult to reproduce and apply the results to specific fracture problems.
To address these problems we present a systematic study of the effects of mechanical layering on the distribution and amount of tensile stress at the tip of a fluid-driven fracture, using the Finite Element Method (FEM). We present a purely static model that considers the conditions for fracture arrest but specifically precludes consideration of the mechanics of propagation. In our models, we consider three variables, namely: (1) Young's modulus contrast (elastic mismatch) across the interface between two layers; that is, the difference in Young's modulus between the layer hosting the fracture tip and the layer ahead of the fracture tip (above the interface); (2) the normalised distance from the fracture tip to the interface (tip proximity to the interface); and (3) the normalised thickness of the layer ahead of the interface. For the latter two parameters, the distance and thickness are normalised to the length of the fracture in the numerical models. One fracture length is 10 m in the models presented. The values are normalised so that the analysis can be applied to fractures of different sizes. We show that all these parameters affect the local tensile stress at the fracture tip, and hence the probability of fracture arrest. Here, the modelling is confined to extension fractures (mode-I fractures). The results therefore have wide implications for the arrest of fluid-driven extension fractures and fracture-related processes such as dyke-fed volcanic eruptions and hydraulic fracturing.
The tensile stress at the tip of the fracture depends on the three variables considered here but also on the dimensions of the fracture. In nature, many mode-I fluid-driven fractures have a length to aperture (thickness) aspect ratio of approximately 1000: 1 ( Gudmundsson, 2011 ;Geshi et al., 2012 ). We therefore use this aspect ratio for the fractures in our models where length refers to the dip dimension of a fracture rather than the strike dimension. Aperture is the size of fracture opening which, for 'frozen' or filled fractures such as dykes and mineral veins, is approximately equal to the measured thickness.

Methods
We used the finite element method software COMSOL Multiphysics to investigate how layering affects the potential arrest of fluid-driven fractures such as dykes and hydraulic fractures. The initial model set up is shown in Fig. 2 . We model the fracture as an elliptical void with an internal overpressure of 0.05 MPa, which is the only loading parameter. The use of overpressure (i.e. pressure above 3 ), rather than total fluid pressure, allows us to take into account the effect of in-situ conditions directly and therefore we do not apply any further loading, either gravitational or tectonic. We used this value of overpressure as it produced tensile stress values at the tip that were in the range of the tensile strengths of many rocks, i.e. 0.5-9 MPa ( Gudmundsson, 2011 ;Amadei and Stephansson, 1997 ). The fracture dimensions were kept constant throughout the study with a fracture length of 10 m and aperture of 0.01 m. These values were chosen as they represent an aspect ratio common for many fluid driven fractures ( Gudmundsson, 2011 ;Geshi et al., 2012 ). The overall model size in vertical cross-section is 300 × 300 m so as to avoid any edge effects on the results, with fixed corners to avoid rigid body translation or rotation. We used a variable mesh with a range of element sizes. The maximum element size is 3 m, remote from the fracture tip, with the minimum element size of 0.006 m being located close to the fracture tip ( Fig. 3 ).
The model fracture is fully confined to what we refer to as the 'host layer', which has the following properties: Young's modulus of 10 GPa, Poisson's ratio of 0.25 and density of 2600 kg m − 3 . These values were chosen as they fall within a realistic range for many sedimentary and igneous rocks ( Gudmundsson, 2011 ) and they also allow for a realistic and large range of Young's modulus values for the layer ahead of the fracture tip to be modelled (i.e. 1.25-80 GPa or a Young's modulus ratio of 0.125-8). Individual layers are mechanically coupled so as to avoid any frictional sliding between the layers. This is a valid assumption for horizontal layers existing at great depths where the normal stress on the interfacial plane is likely to be equal or close to the maximum prin- cipal compressive stress ( Boersma et al., 2020 ). A contrast or change in either Poisson's ratio or density over any realistic range is likely to have a negligible effect on stress magnitudes in comparison to a contrast or a change in Young's modulus, as suggested by earlier studies ( Le Corvec et al., 2018 ;Gudmundsson et al., 2013 ) and supported by our preliminary numerical modelling results. Consequently, we did not conduct an extensive study into the effects of Poisson's ratio or density contrast on tensile stress at the fracture tip.
Unless stated otherwise, the layer ahead is 10 m thick, and the distance between the fracture tip and the interface with the layer ahead is also 10 m. As the fracture length is 10 m, both this distance and the thickness of the layer ahead are equal to 1 fracture length. We consider that the results from this study are likely scale dependant. For this reason, we also quote thicknesses and distances in terms of fracture length. The results can then more easily be applied to hydraulic fractures over a wide range of sizes.
When modelling how the thickness of the layer ahead may affect the tensile stress at the fracture tip we considered a range of thicknesses for that layer, i.e. between 1 m and 80 m, as this adequately covers the range of unit thicknesses found in many sedimentary and volcanic sequences. This also represents a range of 0.1-8 fracture lengths. In terms of hydraulic fractures associated with oil and gas reservoirs, this scale is likely to be a good representation. However, for dyke/sheet intrusions, which can be orders of magnitude larger, a smaller lower bound limit may be more representative to simulate, for example, fine ash layers. Due to computational constraints, we were unable to model thinner layers as they require substantially smaller mesh element sizes.
A range of 0.125-64 m from the fracture tip to the layer ahead was used when considering how a change in proximity to the interface ahead may affect the tensile stress at the fracture tip. This represents a range of 0.0125-6.4 fracture lengths. Again, this constitutes a realistic range for hydraulic fractures propagating in both oil and gas reservoirs and volcanic systems. But for dyke emplacement, a range including a lower limit might be more representative since these are typically much larger. For reasons already discussed, it was not possible to model the fracture tip at such close proximities to the layer ahead due to the same computational constraints.
In the case where the Young's modulus ratio is 1 (i.e. the Young's modulus of the host layer and the layer ahead are the same), the model can be considered homogeneous, and as such any change in the layer ahead (position or thickness) should not affect the tensile stress at the fracture tip. However, we found that when modelling the effect of changing all three parameters simultaneously, there were small differences (mean of < 0.1% difference across all model runs) in the tensile stress in this homogeneous case. These small differences are caused by an automatic fining of the mesh when the fracture tip is closer to the layer ahead, and as the thickness of the layer ahead is reduced. We therefore normalised the results so that, for the homogeneous case, the tensile stress was always the same. This allowed for direct comparisons between model runs.
The absolute values in this study are only valid for an example which matches the model parameters. For this reason, we compare the results using a scale related to the length of the fracture, but also report percentage changes in the tensile stress at the fracture tip.

Results and discussion
Here we present results from 529 model runs where we investigated the three variables, namely Young's modulus contrast, layer thickness, and distance from fracture tip to interface. We analyse the variables first individually and then combined.

Effect of Young's modulus contrast
We modelled how the tensile stress at the tip of a fluid-driven fracture was affected by a succeeding layer (referred to as the layer ahead) with a Young's modulus that was different from or contrasting to the layer hosting the fracture (referred to as the host layer). The properties of the host layer remained constant throughout the analysis, namely a Young's modulus of 10 GPa, a Poisson's ratio of 0.25, and a density of 2600 kg m − 3 . Fig. 4 shows the variation in fracture-tip tensile stress as a function of Young's modulus ratio between the host layer and the layer ahead (A), along with model extracts showing the tensile stress distribution when the Young's modulus of the layer ahead was 1.25 GPa (a ratio of 0.125) (B) and 80 GPa (a ratio of 8) (C). Clearly, the contrast in Young's modulus between the host layer and the layer ahead affects the tensile-stress field at the tip of the hydraulic fracture. When the Young's modulus of the layer ahead and the host layer are the same, i.e. the model is homogeneous, the tensile stress at the fracture tip is 3.84 MPa. When the Young's modulus of the layer ahead is lower than that of the host layer, producing a Young's modulus ratio of < 1, the tensile stress at the fracture tip is higher than in the homogeneous case. This is as expected because a more compliant material with a lower Young's modulus is less able to concentrate stress, and so the tensile stress concentrates in the (stiffer) host layer. Therefore, for lower Young's modulus ratios the tensile stress at the fracture tip is higher. Conversely, when the Young's modulus of the layer ahead is higher than in the host layer, producing a Young's modulus ratio of > 1, the tensile stress at the fracture tip is lower than in the homogeneous case. This is because the tensile stress becomes concentrated in the stiffer layer ahead (see Fig. 4 C). In both Fig. 4 B and C the tensile stress drops significantly within a few tens of centimetres from the fracture tip (i.e. from approximately 4 MPa to less than 0.1 MPa). The variation of tensile stress at the fracture tip with respect to the variation in Young's modulus ratio follows or fits a power-law function ( Fig. 4 A). For the highest Young's modulus ratio modelled (i.e. 8) the tensile stress at the fracture tip is 13% lower than that in the model with the lowest Young's modulus ratio (i.e. 0.125).
The results suggest that, for any given value of host-rock tensile strength, a larger fluid overpressure would be required to propagate a fracture when the layer ahead is stiffer than the host layer. This is in excellent agreement with analytical studies which indicate that many extension fractures become arrested at interfaces where the layer ahead of the interface is stiffer than the layer below the interface that hosts the propagating fracture ( Gudmundsson, 2011 ;Gudmundsson et al., 2010 ;Barnett and Gudmundsson, 2014 ) (i.e. the elastic mismatch mechanism described in the Supplementary Information)

Changes in thickness of the layer ahead
To extend our analysis, we now consider how changes in both the Young's modulus ratio and the thickness of the layer ahead affect tensile stress at the fracture tip. We do this by running two sets of models: (1)

Fig. 5. Variation in the maximum tensile stress
T at the fracture tip with changes in thickness of the layer ahead, for Young's modulus ratios of 0.5 and 2. A power-law function has been fitted to both datasets, and the tensile stress for the homogeneous case is plotted as a dashed red line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) models with a Young's modulus ratio of 0.5 where we vary the thickness of the layer ahead from 2 m to 80 m, and (2) models with a Young's modulus ratio of 2 where we also vary the thickness of the layer ahead from 2 m to 80 m. With both these modulus ratios, the layer ahead has a Young's modulus contrast of 2. As in the previous model runs, the distance between the fracture tip and the interface remains at 10 m (equal to 1 fracture length).

Young's modulus ratio of 0.5
When the Young's modulus of the layer ahead is lower than that of the host layer, the concentration of tensile stress at the fracture tip is higher than that of the homogeneous case. This effect is greater when the more compliant layer ahead is thicker. Thus, the maximum tensile stress at the fracture tip increases as a (positive) power-law function of the increase in thickness of the layer ahead ( Fig. 5 ).

Young's modulus ratio of 2
When the Young's modulus of the layer ahead is higher than that of the host layer, the concentration of tensile stress at the fracture tip is lower than that of the homogeneous case. This is because much the stress becomes concentrated in the stiffer layer ahead. Consequently, the tensile stress at the fracture tip is lower when the comparatively stiff layer ahead is thicker. The maximum tensile stress at the fracture tip decreases as a (negative) power-law function of the increase in thickness of the layer ahead ( Fig. 5 ).
In both implementations of the model (ratios of 0.5 and 2) the tensile stress at the fracture tip approaches the homogeneous value of 3.84 MPa as the thickness of the layer ahead approaches zero ( Fig. 5 ). Thus, as the thickness of the layer ahead is decreased, its effect on the tensile stress at the fracture tip also decreases. When the layer ahead has the maximum modelled thickness of 80 m, the peak tensile stress is ~4% higher than the homogeneous value of 3.84 MPa for a modulus ratio of 0.5, and ~4% lower for a modulus ratio of 2.

Change in distance between the fracture tip and the interface
Here we consider how changes in the distance from the fracture tip to the interface affect the tensile stress concentration at the fracture tip for modulus ratios of 0.5 and 2. In these model runs, the thickness of the layer ahead is kept constant at one fracture length (10 m).

Young's modulus ratio of 0.5
Here the tensile stress at the fracture tip is significantly higher than for the homogeneous case when the fracture tip is close to the interface with the more compliant layer ahead (solid symbols in Fig. 6 ). When the fracture tip is close to the interface, there is less host-layer material between the fracture tip and the interface to dissipate the stress, resulting in a greater tensile stress concentration at the fracture tip. Here, the increase in tensile stress with decreasing distance to the interface is reasonably well fit with a power law.

Young's modulus ratio of 2
Under these conditions, the tensile stress at the fracture tip is significantly lower than for the homogeneous case when the fracture tip is close to the interface with the stiffer layer ahead (open symbols in Fig. 6 ). When the fracture tip is closer to the interface with the stiffer layer, it concentrates more of the stress, thereby reducing the stress at the fracture tip. Here, the decrease in tensile stress with decreasing distance to the interface is also reasonably well fit with a power law.
For both conditions, the model demonstrates that the tensile stress at the fracture tip is highly dependant on the proximity of the fracture to the layer ahead. However, this dependence decreases significantly as the distance between the fracture tip and the layer ahead is increased. For both conditions, our data are reasonably well fit by a power law up to the distance where the tensile stress at the fracture tip is the same value as the homogeneous case. At greater distances, it is a poor fit as our values for both ratios are the same as the homogeneous case. The results show that there is essentially no effect on fracture-tip tensile stress in Fig. 6. Variation of the maximum tensile stress T at the fracture tip with distance from the interface, for Young's modulus ratios of 0.5 and 2. A powerlaw function has been fitted to both sets of model output data, and the tensile stress for the homogeneous case is plotted as a dashed red line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the far-field, where the distance to the interface exceeds approximately 30 m, or 3 fracture lengths Furthermore, the magnitude of change in fracture-tip stress depends on the modulus contrast. For a fracture close (0.125 m) to an interface with a more compliant layer, where the modulus ratio is 0.5, the fracture-tip tensile stress is 37% higher than in the far-field. By contrast, for a fracture approaching close to an interface with a stiffer layer, with a modulus ratio of 2, the fracture-tip tensile stress is 26% lower than in the far-field.
In a geological context, these results suggest that, for any given loading condition, when the tip of a fluid driven fracture approaches a layer interface, the tensile stress at the fracture tip gradually increases if the layer being approached is more compliant than the host layer, but gradually decreases if the layer being approached is stiffer than the host layer. Hence, such fractures (for example, dykes) propagating through a comparatively compliant layer but approaching the interface with a comparatively stiff layer may become arrested well before they reach the interface -a phenomenon commonly observed in the field ( Gudmundsson, 2011 ;Barnett and Gudmundsson, 2014 ) ( Fig. 1 ). By contrast, a fracture propagating through a comparatively stiff layer but approaching the interface with a comparatively compliant layer is likely to become arrested only at the actual interface. This interpretation is, again, well supported by field observations which show, for example, that most fractures in stiff limestone layers propagate right up to the interfaces with adjacent and more compliant shale layers ( Gudmundsson, 2011 ).
Similar observations have been reported from uniaxial compression tests on samples with multiple layers of granite, sandstone and siltstone ( Douma et al., 2019 ). In these experiments, the samples comprised three rock layers, with the top and bottom layer being of the same material and the middle layer being a contrasting material. Douma et al. (2019 ) characterised the contrast in mechanical properties using the unconfined compressive strength (UCS) of the different materials, and concluded that when the contrast was greatest between the layers fractures propagating from the weaker layer would sometimes become arrested before reaching the contact with the stronger layer. Although their work was based on a comparison of the UCS values of the contrasting layers rather than Young's modulus, we consider these results to be a useful comparison with our work since the two parameters are closely related ( Douma et al., 2019 ). Furthermore, the authors also suggest that as the fractures in their experiments approach the interface the stress intensity factor at the crack tip approaches zero. As the mode-I stress intensity factor and tensile stress at the crack tip are also closely related ( Forbes Inskip et al., 2018 ;Gudmundsson, 2011 ;Chandler et al., 2016 ;Zhang, 2002 ), their results generally support our observations and interpretations.

Varying all three parameters simultaneously
Of the three parameters discussed here that influence the fracturetip tensile stress, the proximity to the interface between mechanically dissimilar layers appears to have the largest effect, while the thickness of the layer being approached appears to have the least effect. Nevertheless, it is important to understand how the tensile stress may vary with combined changes in all three parameters In order to combine changes in all three parameters simultaneously, we ran a series of multi-parametric model runs with (1) Young's modulus ratios between 0.0125 and 8, (2) fracture tip to interface proximities between 0.125 m and 40 m (0.0125 to 4 fracture lengths), and (3) layer thicknesses of 1 m and 80 m (0.1 and 8 fracture lengths). The outputs from all model runs are provided in the Supplementary Information. A representative selection of the outputs is displayed in Fig. 7 .
When combining the variations in all three parameters, we find that the tensile stress at the fracture tip can range from a minimum of 1.81 MPa to a maximum of 11.67 MPa; or approximately 0.5 to 3 times that for the homogeneous case ( Fig. 6 ). Unsurprisingly, this range of values is significantly higher than when each of the variables is changed independently. Both the maximum value (approximately 3 times higher than for the homogeneous case) and the minimum value (approximately half that for the homogeneous case) occur when the fracture tip is clos- Fig. 7. Tensile stress at the fracture tip as a function of variations in all three parameters: (1) Young's modulus ratios between 0.0125 to 8), (2) interface proximities from 0.125 m to 40 m, and (3) succeeding layer thicknesses of 1 m and 80 m. The tensile stress for the homogeneous case is plotted as a horizontal dotted black line. All other lines plotted are power-law fits for their respective data set.
est to the interface and when the thickness of the layer ahead is greatest (80 m or 8 fracture lengths).
The highest and lowest tensile stresses occur at the model parameter extremes, and while the lower values of Young's moduli used in this study are common in terms of up-scaled rock masses, some of the higher values (i.e. > 20 GPa) may be uncommon, or even unrealistic, for up-scaled fractured rock masses. A comprehensive study by Heap et al. (2020 ) reviewed existing Young's moduli data for volcanic rocks, and showed that even at depth where Young's moduli values are likely to be higher than at the surface, volcanic rocks rarely exhibit intact Young's moduli greater than 30 GPa. Furthermore, when upscaling these values to take into consideration the presence of macroscopic natural fractures, Heap et al. (2020) found that realistic Young's moduli values for volcanic rocks were even lower and rarely above 10 GPa. The intact Young's moduli of some sedimentary rocks (in particularly carbonates) and metamorphic rocks has been shown to be as high as 80 -100 GPa ( King, 1983 ;Eissa and Kazi, 1988 ;Brotons et al., 2015 ). However, these are Young's moduli values of intact rocks (i.e. taken from plug scale intact rock samples) and by extending the work of Heap et al. (2020) to other rock types, up-scaled values are likely to be somewhat lower. Bearing this in mind, we nevertheless stress that it is the contrast in Young's modulus between the individual layers (rather than the absolute moduli values) that is the key parameter to consider in our analysis. Therefore, even though our study contains model iterations with Young's moduli values of 10 GPa and 80 GPa for the host layer and layer ahead, respectively, which may be uncommon for a volcanic complex in situ, a Young's modulus ratio of 8 is not uncommon.
For model iterations with less severe variations in Young's modulus ratio, layer thickness and interfacial proximity, the variations in tensile stress at the fracture tip are smaller but remain significant. For example, when the thickness of the layer ahead is only 1 m, and the fracture tip is 0.125 m from the interface, we still observe maximum and minimum tensile stresses that are 66% higher (6.36 MPa) and 41% lower (2.28 MPa) than the homogeneous value, and thus still significant. For example, in the Vaca Muerta formation of Argentina (a layered uncon-ventional hydrocarbon reservoir, where layers of organic-rich shale are interbedded with other units such as limestone, ash beds and sills), individual layer thicknesses vary from a few centimetres to a few metres ( Sosa et al., 2017 ). Therefore, if a fluid-driven fracture of length < 10 m were modelled assuming a homogeneous formation with no mechanical contrast between layers, the resulting stress magnitudes would likely contain significant errors when compared with a more realistic, heterogeneous layered model.
The same is true for crustal segments hosting volcanoes where individual (commonly compliant) pyroclastic layers can range in thickness from a few tens of centimetres to tens or even hundreds of metres ( Gudmundsson, 2020 ), whereas (significantly stiffer) sills and lava flows may be as thick as tens or (particularly sills) hundreds of metres. The tensile stress at the tip of a dyke with a height (dip dimension) of < 1 km may be very different depending on whether the hosting crustal segment is modelled as being homogeneous or heterogeneous (i.e. containing layering). Stratovolcanoes normally contain many compliant pyroclastic (tuff) layers alternating with much stiffer lava flows (and some sills) through which dykes must propagate in order to reach the surface to erupt ( Gudmundsson, 2009 ). Most dykes, however, do not reach the surface but rather become arrested at some depth in the volcano ( Gudmundsson, 2011 ;Barnett and Gudmundsson, 2014 ;Geshi et al., 2012 ;Gudmundsson and Philipp, 2006 ). Based on the model presented here, dyke arrest when the layer ahead (above the interface) is more compliant than the layer hosting the dyke tip is most likely to occur directly at the interface. This follows because the tensile stress at the fracture tip increases as it approaches the interface, encouraging further propagation. This comparatively high tensile stress may also help to open the interface, resulting in a Cook-Gordon delamination ( Gudmundsson, 2011 ;Barnett and Gudmundsson, 2014 ). By contrast, dyke arrest when the layer ahead is stiffer than the host layer can occur either at the interface itself, or at some distance before the interface is reached. This occurs because the tensile stress at the dyke tip decreases as it approaches the interface. Thus, the arrest location is not confined to the interface itself, but actually becomes increasingly more likely as Here the fracture may become arrested before reaching the interface between the two layers. (B) Fracture is propagating from a relatively stiff layer towards a comparatively compliant layer, where the fracture length is approximately the same as the thickness of the layer ahead. Here the fracture is more likely to reach, and be arrested at, the interface between the two layers. (C) Fracture has originated from a relatively compliant layer compared to the comparatively stiff layer above, but where the fracture length is much greater than the thickness of the layer ahead. Here the effects of the stiffer layer ahead in reducing the tensile stress at the fracture tip are overcome by the length of the fracture being much greater than the thickness of the layer ahead and, as such, the fracture is more likely to propagate from the compliant layer into the stiff layer.
the dyke tip approaches the interface. This is, indeed, what is commonly observed in the field ( Fig. 1 ). Furthermore, the reduction in the fracture tip tensile stress as the tip comes into closer proximity with the interface with a stiffer layer is one of the primary reasons why elastic mismatch commonly leads to fracture arrest. Thus, the present numerical model provides strong support for the elastic mismatch mechanism of fracture arrest  which is based on analytical considerations, but only considers arrest at an interface. More specifically, our results indicate that when mode I (extension) fractures approach an interface, arrest of the fracture tip is most likely when the layer ahead of the interface is stiffer than the layer hosting the fracture tip ( Fig. 8 ).
The results of this study demonstrate the importance of correctly identifying layering within geological sequences in order to understand the likelihood of fracture arrest. Specifically, it is important to ascertain individual layer thicknesses, the Young's moduli of the individual layers and the proximity of propagating fractures to the layers ahead (i.e. the interfacial distance). In reality it can be difficult to determine these parameters. However, using a combination of field and laboratory methods it is possible to determine these parameters to at least a first order approximation. For example, the Young's modulus of the individual layers (and therefore the modulus contrast between layers) can be determined using standard laboratory experiments on samples gathered from core or outcrop material ( Heap et al., 2020 ;Heap and Faulkner, 2008 ;ISRM, 1970 ). The dynamic Young's modulus can also be calculated from velocity data gathered using sonic well log tools ( Heap et al., 2020 ;Eissa and Kazi, 1988 ;Brotons et al., 2015 ). As it is the Young's modulus contrast between layers that is important rather than Young's modulus of individual layers, as a first order approximation either dynamic or static Young's moduli values can be used (as long as static values are compared to static values and dynamic values are compared to dynamic values). Fracture locations can be determined using microseismic monitoring in actively producing hydrocarbon reservoirs ( Maxwell and Norton, 2012 ;Majer et al., 2007 ;Karamzadeh et al., 2019 ). Although this may not give a precise location of the fracture tip, continuous monitoring during stimulation treatments will provide an approximate location. Individual layer thicknesses can be determined from either core or well log data, should they exist, or from seismic survey data. Unfortunately it can be difficult to distinguish individual layers from seismic survey data alone, unless the layers are particularly thick or have a significant velocity contrast to adjacent layers ( Widess, 1973 ;Mavko et al., 2009 ;Chung and Lawton, 1999 ). Muon tomography can also be used to understand the structure of layered sequences, and has been used to map the structure of active stratavolcanoes ( Nishiyama et al., 2014 ;Lesparre et al., 2012 ;Le Gonidec et al., 2019 ).

Conclusions
This study presents outputs from a series of numerical models which highlight how crustal layering affects fracture propagation. Our models simulate fluid-pressure driven mode-I fractures which are analogous to volcanic intrusions such as dykes and man-made hydraulic fractures which are often generated to stimulate oil and gas reservoirs. The results, however, may be applicable to other modes of fracture propagation, with some modifications.
In our models, we altered the thickness of the layer ahead of a propagating fracture, the proximity of the fracture tip to the interface with the layer ahead, and the Young's modulus contrast between the layer hosting the fracture and the layer ahead. The results show that all three parameters affect the tensile stress concentration at the fracture tip, and therefore also influence the potential for fracture propagation or, alternatively, arrest. Overall, the model results show that: • Tensile stress at the fracture tip increases as the tip approaches an interface if the layer ahead (above the interface) is more compliant than the layer hosting the fracture tip, but decreases when the layer ahead is stiffer than the hosting layer. • The contrast in Young's modulus between the layers (elastic mismatch) has a larger effect on the fracture-tip tensile stress as the layer ahead of the propagating fracture increases in thickness and as the proximity of the fracture tip to the interface decreases. • Fracture arrest is most likely to occur exactly at the interface when the layer ahead is more compliant than the host layer. By contrast, when the layer ahead is stiffer than the host layer fracture arrest can occur either at the interface, or at some distance within the host layer before the interface is reached. These results are in agreement with field observations of arrested extension fractures.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.