Fracturing and Dome‐Shaped Surface Displacements Above Laccolith Intrusions: Insights From Discrete Element Method Modeling

Inflation of viscous magma intrusions in Earth's shallow crust often induces strain and fracturing within heterogeneous host rocks and dome‐shaped ground deformation. Most geodetic models nevertheless consider homogeneous, isotropic, and linear‐elastic media wherein stress patterns indicate the potential for failure, but without simulating actual fracturing. We present a two‐dimensional Discrete Element Method (DEM) application to simulate magma recharge in a pre‐existing laccolith intrusion. In DEM models, fractures can propagate during simulations. We systematically investigate the effect of the host rock toughness (resistance to fracturing), stiffness (resistance to deformation) and the intrusion depth, on intrusion‐induced stress, strain, displacement and spatial fracture distribution. Our results show that the spatial fracture distribution varies between two end‐members: (a) for high stiffness or low toughness host rock or a shallow intrusion: extensive cracking, multiple vertical surface fractures propagating downward and two inward‐dipping highly cracked shear zones that connect the intrusion tip with the surface; and (b) for low stiffness or high toughness host rock or a deeper intrusion: limited cracking, one central vertical fracture initiated at the surface, and two inward‐dipping fractures at the intrusion tips. Abrupt increases in surface displacement magnitude occur in response to fracturing, even at constant magma injection rates. Our modeling application provides a novel approach to considering host rock mechanical strength and fracturing during viscous magma intrusion and associated dome‐shaped ground deformation, with important implications for interpreting geodetic signals at active volcanoes and the exploitation of geothermal reservoirs and mineral deposits.

Owing to the lack of direct observations at active volcanic systems, the factors controlling the emplacement and inflation of thick sills and laccoliths have been instead studied at exhumed, inactive volcanic plumbing systems.Large strains have been found to lead to pervasive elastic bending, inelastic compaction, and disruption by fracturing and faulting, resulting in so-called "forced folds" and sequences of normal faults (Bischoff et al., 2019;Burchardt et al., 2019;Magee et al., 2018;Mattsson et al., 2018Mattsson et al., , 2020;;Reeves et al., 2018;P. I. Wilson et al., 2016).Importantly, tilted eruptive cones and overburden sequences, brittle deformation within the viscous magma, and variations in the geochemistry of the magma indicate that a rapid initial phase of magma-driven fracture opening is often followed by inflation into a thick intrusion (Burchardt et al., 2019;Mattsson et al., 2018;Morgan et al., 2008;Petronis et al., 2019).The result is often a complex pattern of fractures and faults in the overlying host rock.
Geological observations from inactive volcanic systems only allow static observations of the post-intrusive state.Scaled laboratory experiments have been used to study the damage mechanisms around the magma propagation front, the transition from a layer-discordant dyke to a layer-parallel sill, and the controls of imposed stress fields on the pathways magma takes to the surface (e.g., Galland et al., 2018;Kavanagh et al., 2018;Menand, 2011;Poppe et al., 2019;Rivalta et al., 2015).Scaled laboratory experiments have also documented how magmatic pressure source depth controls the geometry of faults formed in granular host materials (Montanari et al., 2017;Warsitzka et al., 2022).Most laboratory simulations encompass a scaled laccolith width of 50-2,000 m, whereas the simulated range of mechanical properties of crustal rock strengths is limited by the available analog granular materials of which the elastic properties remain unconstrained: cohesion values between 10 5 and 10 7 Pa, tensile strengths between 10 4 and 10 7 Pa, and friction angles between 30 and 40°(e.g., Galland et al., 2018;Merle, 2015;Poppe et al., 2021).Stresses cannot be quantitatively measured.
Analytical and numerical models are more versatile in investigating the dynamic process of magma intrusion.These models have focused on the effect of topographic loads, host rock properties, and magma rheology on the deformation of the overlying host rocks.Most analytical and Boundary Element Method (BEM) models consider a linear-elastic response of the shallow crustal rocks to the intrusion pressure, and therefore indicate the potential for failure in the host rock (Bunger & Cruden, 2011;Galland & Scheibert, 2013;Grosfils et al., 2015;Michaut, 2011;Morgan, 2018).The Finite Element Method (FEM) considers strain and stress patterns to determine fracture pathways (Rivalta et al., 2019).Although these studies have provided essential insights into how the inflation and position of magmatic pressure sources in the subsurface affect ground deformation patterns over time, the assumption of linear elasticity prevents consideration of the effect of fracturing on the stress field.This also limits the exploration of progressive weakening that can be expected from the accumulation of smallscale cracking around larger-scale fractures.Recent models that include an inelastic component, have shown that inelastic deformation, host rock heterogeneities, and far-field tectonic stresses can modify how magma intrusions grow and displace the surface (e.g., Carrier et al., 2015;Gerbault et al., 2018;Gill & Walker, 2020;Got et al., 2019;Haug et al., 2017;Souche et al., 2019).These models dominantly focus on the initial emplacement of viscous intrusions, whereas existing field observations show that later-stage inflation may affect the eventual failure of the overlying host rocks (Burchardt et al., 2019;Mattsson et al., 2020).However, the required remeshing of often large elements to simulate large-strain problems in FEM is computationally prohibitive for crustal-scale models.
The Discrete Element Method (DEM) is an alternative numerical method, whereby host rock is represented as spherical particles that interact at inter-particle contacts.Breakage of these contacts therefore represents cracking of the host rock, allowing DEM models to explore damage accumulation and fracture propagation during addition of material to the system, that is, magma intrusion.DEM models have been used to systematically investigate the effects of variable host rock mechanical properties during volcanic dome growth (Harnett et al., 2018) and instability (Harnett et al., 2022), depletion-induced collapse at sinkholes (e.g., Al-Halbouni et al., 2018;Hafver et al., 2014;Smart et al., 2011) and calderas (e.g., Harnett et al., 2023;Holohan et al., 2017;Woodell et al., 2023), and dyke-fed sandstone intrusions (Meng & Hodgetts, 2020).Zhao et al. (2008) have used DEM models to investigate fractures above laccolith intrusions in homogeneous and 2-layered models, but they did not study the effect of mechanical host rock properties on the observed fracture patterns and displacements.
We present results from a 2D DEM application of magma influx into a laccolith-shaped intrusion that allows the systematic study of the effect of host rock mechanical properties and intrusion depth on dynamic fracturing and surface displacement patterns.We systematically vary the numerical parameters and relate those to a realistic range of mechanical rock properties to independently vary the host rock ability to resist bending (stiffness) and fracturing (toughness).Our study benefits of the unique feature of DEM models wherein we can track the effect of the dynamic fracturing of the host rock on principal stresses, normalized finite shear strain, and surface displacements.

The Discrete Element Method
The Discrete Element Method (DEM) was first developed to model granular media (Cundall & Strack, 1979).DEM discretizes the modeled material into independent rigid particles with a finite mass.In a two-dimensional (2D) model, these rigid particles are disk-shaped.The particles can move independently according to Newton's laws of motion and interact with other particles at contact points.Per model time step, a force and moment resulting from both gravitational and contact forces is calculated per particle to obtain their new positions and velocities.We use the DEM Particle Flow Code 2D software (PFC2D) from Itasca Consulting Group, Ltd. to implement our model application (https://www.itascacg.com/software/PFC).
Properties of contacts between particles can be defined so that the discontinuous assemblage of particles acts as a continuous block of rock that can resist tension and shear stresses (Potyondy & Cundall, 2004).DEM particles should not be seen as individual rock grains or blocks, but rather as discrete packets for the purpose of computation.By using contact models provided by PFC, all mechanical properties of the modeled material can be controlled, such as the failure envelope, elastic parameters, and the inter-particle friction.The principal model parameters are summarized in Table 1.The used contact models are defined in such a way that contacts are considered, and solved, along virtual continuous lines equal to the smaller diameter of the two particles in contact (Itasca, 2021;Potyondy & Cundall, 2004).This contact definition also ensures that the developed stresses are independent of particle size.Particles are bonded together until a predefined bond tensile strength or bond shear strength (controlled by the bond cohesion) is reached.Unbonded particles (i.e., post-failure) then mainly interact with one another through friction and cannot resist tensile forces.This way, failed bonds, or "cracks," can align to form fractures that can propagate dynamically throughout the medium as it deforms.

The Magma Intrusion Model
We simulate the injection of new magma into a pre-existing shallow intrusion and explore the effect on the deformation of the host rock as a function of host rock properties and intrusion depth.Overall, laccolith intrusions are found in the shallowest 3 km of the Earth's crust (Cruden et al., 2018).We vary intrusion depth in our models between 500 and 1,000-1,500 m similar to proposed depths for the tops of well-studied laccoliths, that is, ≈500 m at Sandfell in Iceland (Mattsson et al., 2018), 500-1,500 m at Kilian and Petit Puy de Dôme, Chaîne des Puys, France (Petronis et al., 2019), and 200-1,500 m of forced fold amplitudes in seismic profiles of the Bight Basin Igneous Complex, Australia (Jackson et al., 2013).
We focus on geometries typical for such rather small-scale laccoliths, including the most studied ones of the Henry Mountains in Utah, USA, which are typically 0.5-2 km wide and 50-500 m thick (Cruden et al., 2018;McCaffrey & Petford, 1997;Morgan et al., 2008).We thus start our simulation with a pre-existing source that is 1,000 m wide and 50 m thick with a flat base and half-elliptical roof shape.The source of deformation is embedded in a 20 km-wide rectangular domain to prevent lateral boundary effects (Figure 1).We use fixed boundary conditions at the model edges to simulate infinite extension of rock laterally and downwards.There is a free surface at the top of the boundary to allow exploration of surface deformation.We place the laccolith at the bottom of the model space and neglect any deformation below the intrusion.The DEM particles have an average radius of 6.65 ± 1.5 m, which is similar or even smaller than ranges used for similar kilometric scale model setups, for example, 5-10 m in Holohan et al. (2015Holohan et al. ( , 2017)), 8-13 m in Woodell et al. (2023), 10-20 m in Smart et al. (2011), 15-25 m in Meng and Hodgetts (2020) or 28.78-43.1 m in Zhao et al. (2008).A particle size ratio of 1.66 is widely used because it ensures that hexagonal packing is avoided and that the particle assemblage can pack to a suitable packing porosity (Al-Halbouni et al., 2018;Holohan et al., 2015Holohan et al., , 2017;;Meng & Hodgetts, 2020;Woodell et al., 2023).Smaller particle radii are not practical due to the high computational cost, which increases non-linearly with a further linear decrease in particle radius.All particles have the same density of 2,700 kg•m 3 .After initial packing, gravity is introduced and a settling routine is performed to optimize the particle packing.Thus the effect of lithostatic pressure is taken into account even if all contacts in the host rock have the same numerical proprieties.Finally, each particle is attributed to one of two groups: rock or magma (respectively gray and pink-red particles in Figure 1).Among the contact models offered in PFC, three different contact models were used to represent as accurately as possible the rock, the magma, and surfaces along produced cracks.Contacts between two host rock particles and between host rock particles and the model walls are defined using the soft-bond contact model (Itasca, 2021;Jiang et al., 2015;Ma & Huang, 2018a, 2018b).In general, particles are allowed to slightly overlap to model pre-failure elastic deformation in an assembly of particles (Cundall & Strack, 1979), as is the case for natural rocks (Jaeger et al., 2009).A soft-bond contact is mainly characterized by a bond effective modulus, a normal to shear stiffness ratio, a bond cohesion, a bond tensile strength, a bond friction coefficient, and softening parameters (Table 1).The other parameters which do not control the particle packing toughness or stiffness follow the values listed in Table 1 or can be inferred from the parameters in it (Itasca, 2021).The main asset of the soft-bond contact is its ability to achieve a realistic value of the ratio of compressive to tensile strength for the simulated rock (Ma & Huang, 2018a, 2018b).When the softening is activated, the bond is not immediately deleted when the bond tensile strength is exceeded.Instead, the bond tensile stress decreases linearly with continued extension according to the softening factor γ until it reaches a threshold percentage of the bond tensile strength, defined by the softening tensile strength factor ζ.This mechanism accounts for the non-linear, inelastic part of the stress-strain behavior observed in natural rock samples as they approach macro-scale failure in laboratory strength tests.We set softening parameters γ = 13.0 and ζ = 0.0 and keep these parameters constant in all simulations to maintain a realistic compressive to tensile strength ratio of approximately 10, as observed for natural rocks (Heap & Violay, 2021).All other main contact parameter values are listed in Table 1.
We approach the magma as an incompressible fluid.Contacts between magma particles are defined with the linear-parallel bond contact model (Potyondy & Cundall, 2004), with no friction, cohesion or tensile strength, and a bond effective modulus of 25 GPa.Its low resistance to relative particle rotation makes this contact model appropriate for simulating fluid-like behavior (Harnett et al., 2018).The values of the other parameters are given in Table 1.The contacts between magma and host rock particles are also defined using the linear-parallel bond contact model, but with a friction coefficient of 0.5.Because the diameter of the particles in our DEM model setup is on the meter-scale, sharp intrusion tips cannot be simulated; our application instead simulates the intrusion of viscous magma with meter-scale propagation fronts, as generally observed in shallow laccolith intrusions.A more precise modeling of the magma fluid dynamics is currently not possible in PFC and will require future work.
When a contact breaks between two host rock particles or between a rock particle and a wall, the soft-bond contact model is replaced with the rolling resistance contact model (Itasca, 2021;Jiang et al., 2015).This contact model is also applied if two rock particles come into contact for the first time such as, for example, following the deformation of the host rock.Parameter values, such as friction coefficient, are the same as those defined for the host rock (Table 1).The main advantage of the rolling resistance contact model is that even if the particles are not bonded, and only interact through friction, the contact resists the relative rotation of particles.This is a realistic representation of the rigid surface contact between two detached blocks.We increase the volume of the intrusion by injecting new magma particles (red particles in Figure 1b) from a conduit at the base of the pre-existing laccolith-shaped intrusion (pink particles in Figure 1).Particles are injected at a constant vertically-upwards rate of 1 m•s 1 .When a particle has completely crossed the permeable portion of the bottom wall, the imposed velocity is removed.Then, the particle velocity is determined only by the forces and moments that result from gravity and particle-particle interactions.The magma intrusion thus is the deformation source and sits within a host rock consisting of rock particles (gray particles in Figure 1).The simulation stops at 100% of injection, which means that the amount of particles in the intrusion is then twice that of the pre-existing intrusion.Rather than explicitly considering time in our model, we consider the proportion of magma injection.Thus, we express the simulation progress in terms of percentage of magma injected and focus on 5%, 25%, and 100% of injection, where 100% represents a doubling of the injected area in the 2D model.Assuming an axial symmetry, the simulation thus represents the doubling of the volume of the initial magma body with ≈0.1 km 3 to a total of ≈0.2 km 3 .Thus 5% and 25% of injection correspond, respectively, to ≈5 × 10 3 km 3 and ≈2.5 × 10 2 km 3 of magma.The total simulated emplacement duration can be approximated to a few days or weeks.

Model Exploitation
Our aim is not to reproduce a particular rock, but to systematically explore the broader numerical parameter space.
Ranges of numerical parameters are given in Table 1 and discussed below.Numerical parameters that control particle bonds are not equal to the bulk mechanical properties of the particle assemblage.We infer the bulk mechanical properties of the modeled rock by performing numerical Uniaxial Compressive Strength (UCS) and Tensile Strength (TS) tests further described in Text S1 in Supporting Information S1 (Huang et al., 2013;Ma & Huang, 2018a).The obtained bulk properties are provided for reference in Table 2 at the end of the results section and can be directly related to mechanical properties of natural rocks.In line with Ma and Huang (2018b); Ma and Huang (2018a), we inhibit shear failure in the rock by setting the bond cohesion to ten times the bond tensile strength.Even if each individual bond breaks in tension, the resulting alignment of cracks can develop shear bands and fractures at the assemblage scale (Huang et al., 2013;Ma & Huang, 2018a).This is illustrated by the numerical UCS tests in which all bonds break in tension, but cracks align along a typical shear fracture (Figure S1 in Supporting Information S1).We focus our study on three key features: rock toughness, rock stiffness, and source depth: 1. Variation in toughness, that is, the rock's resistance to fracturing, is mainly controlled by the bond tensile strength and the bond cohesion (at a 1:10 ratio).For a higher toughness, a rock can accommodate more elastic deformation before breaking.2. Variation in rock stiffness, that is, the rock's resistance to elastic deformation, is mainly controlled by the bond's effective modulus and normal to shear stiffness ratio.The higher the stiffness of the rock, the more force is required to obtain a given amount of deformation.3. Variation in source depth has been shown to affect surface displacement patterns for simple pressure sources in a homogeneous linear-elastic half-space (Battaglia et al., 2013;Masterlark, 2007).

Results Processing and Visualization
For each particle, the normalized finite shear strain is determined using the Cauchy-Green deformation tensor determined from the displacement gradient tensor (Cardozo & Allmendinger, 2009;Schöpfer et al., 2006).The complete method is described by Harnett et al. (2018).Due to large particle displacements near cracks, some particles have an anomalously high normalized finite shear strain compared to most other particles.However, considering all simulations at the same time, at least 90% of the finite shear strain data are between 0.0% and 0.05%.To help visual comparisons between different simulations, we apply the same cutoff criterion of 0.05% to all plots and subplots.All finite shear strains above a value of 0.05% are set to 0.05%.Finite shear strains are then normalized to emphasize the relative shear strain distribution within one simulation.
We visualize breakage of individual bonds (i.e., cracks) by a black line perpendicular to the broken bond.These lines have the same length as the corresponding deleted contact.We discuss damage accumulation by considering the amount of broken contacts for a given simulation state, expressed as a percentage of initially intact contacts.
Similarly to the normalized finite shear strain, principal stresses are determined for each particle.The average stress tensor of each individual particle is diagonalized to determine the particle's principal stress magnitudes as in Holohan et al. (2015).Compressive stress is negative and tensile stress is positive.In the highly cracked zones in the host rock, some totally detached particles or small particle chains locally concentrate very high principal stress values compared to others.To better visualize and compare the principal stress patterns, a same cutoff was applied  To allow comparison between different simulations at each percentage of injection, we normalize the vertical displacements by the maximum vertical displacement detected at 100% of injection of the corresponding simulation.To determine the lateral extent of the vertically displaced surface, we include all surface particles of which their absolute vertical displacement is higher than 0.5% of the value of the maximum vertical displacement detected in the time step at 100% of the corresponding simulation.

Results
In this section, we present the effect of the host rock toughness and stiffness and the source depth on finite shear strain and principal stresses.For each performed simulation plots of the raw model output in terms of particle position at 5%, 25% and 100% are given in Figures S3 to S15 in Supporting Information S1.

Effect of Rock Toughness and Stiffness
We first vary the toughness.For selected representative simulations, the normalized finite shear strain and crack distribution are shown for different host rock toughnesses at 5%, 25% and 100% of injected 2D area (Figure 2).We modeled bond tensile strengths of 1 MPa (Figures 2a-2c), 5 MPa (Figures 2d,-2f), and 15 MPa (Figures 2g-2i).The bond cohesion is ten times larger than that of the bond tensile strength in all simulations, as explained in the Method Section above.Movie files of the temporal evolution of damage and strain along with horizontal and vertical surface displacement for all our simulations are provided in Supporting Information S1.The simulation with a bond tensile strength of 1 MPa represents an intrusion into a host rock with a low toughness (Figures 2a-2c).At 5% of injection, the host rock is already cracked and a zone of high normalized finite shear strain has developed near the intrusion with its maximum above the intrusion center (Figure 2a).As the injection progresses to 25%, host rock cracking increases, the zone of high normalized finite shear strain grows outward and upward into two inclined zones that reach the surface, and several tensile fractures open at the surface (Figure 2b).At 100% of injection, the host rock is highly fractured, more tensile fractures have opened at the surface, and there is increased normalized finite shear strain accumulation in the intact host rock (Figure 2c).
Simulation results for a host rock with an intermediate toughness and a bond tensile strength of 5 MPa are shown in Figures 2d-2f.By 5% of injection, cracks have opened at the center of the intrusion roof and started to define a short vertical fracture zone propagating toward the surface (Figure 2d).Then, by 25% of the total injection volume, a nearly vertical tensile fracture has opened at the surface above the center of the intrusion, and has propagated downward and joined the upward propagating fracture zone (Figure 2e).Finally, at 100% of injection, the fracture pattern has become asymmetric, with cracks accumulated at the intrusion edges and forming two fractures that bound a zone of high normalized finite shear strain on one side of the vertical fracture (Figure 2f).The aperture of the surface fracture has increased from 13 to 46 m between 25% and 100% of injection.
A bond tensile strength of 15 MPa is representative of a host rock with a high toughness (Figure 2g-2i).The host rock is barely fractured at 5% of the injection, and a zone of high normalized finite shear strain has developed near the intrusion with its maximum above the intrusion center (Figure 2g).By 25% of injection, the host rock remains mostly intact.The extent of the zone of high normalized finite shear strain grows centrally above the intrusion in one vertical elliptical zone, and one small zone of high normalized finite shear strain is visible at the surface (Figure 2h).At 100% of injection, the fracture pattern is concentrated into a single 48 m-wide central fracture which opens from the surface downward, and two less developed inward-dipping inclined fractures that propagate from the intrusion edges (Figure 2i).
Generally, Figure 2 shows that higher host rock toughness requires a greater injected volume to open tensile fractures at the surface (Figures 2a, 2e, and 2i).Higher toughness results in less extensive fracturing of the host rock and more localized deformation overall (Figures 2c, 2f, and 2i).
Regardless of the host rock toughness, when vertical fractures open at the surface, we find a reduction in the surrounding normalized finite shear strain (Figures 2a,2b,2c,2e,2f,and 2i) and the high normalized finite shear strain concentrates in two zones that incline away from the intrusion edges (f and i in Figure 2).
Next, we vary the host rock stiffness, that is, its ability to resist elastic deformation.Figure 3 shows the normalized finite shear strain and fracture patterns for selected representative simulations for different host rock stiffnesses, for 5%, 25%, and 100% of magma injection.The bond effective modulus of the selected simulations was varied between 3 GPa (Figures 3a-3c), 15 GPa (Figures 3d-3f), and 30 GPa (Figures 3g-3i).The bond tensile strength was kept constant at 5 MPa and the bond cohesion at 50 MPa.
The simulation with a bond effective modulus of 3 GPa and a bond tensile strength of 5 MPa is the host rock with the lowest stiffness (Figures 3a-3c), and is the same simulation as that used for the intermediate toughness rock described above (Figures 2d-2f).
The simulations with an intermediate stiffness (effective modulus = 15 GPa) (Figures 3d-3f) and a high stiffness (effective modulus = 30 GPa) (Figures 3g-3i), show similar spatial distribution of cracks along with patterns and magnitudes of normalized finite shear strain.A fractured zone develops vertically above the intrusion center in the initial 5% of injection (Figures 3d and 3g).Fractures have appeared on each side of the vertical fractured zone by 25% of injection (Figures 3e and 3h).Finally, by 100% of injection, the host rock is highly fractured and multiple tensile fractures have opened at the surface (Figures 3f and 3i).Two inward-dipping fractures curve at 200 m below the surface and clearly developed fractured zones connect the intrusion edges and the surface (Figures 3f  and 3i).In stiffer host rock (effective modulus = 30 GPa), the two inclined fractured zones and surface fractures develop earlier and are visible at 25% of injection (Figure 3h).
Overall, in stiffer host rocks, tensile fractures open at the surface at earlier stages of injection (Figures 3b,3d,3g,and 3h) and the host rock is more extensively fractured (Figures 3c,3f,and 3i).At 100% of injection, stiffer host rock results in both a lower magnitude and narrower distribution of normalized finite shear strain in the rock (Figures 3c, 3f, and 3i).This is because the strain is accommodated by increased fracturing.

10.1029/2023JB027423
For the same simulations as in Figures 2 and 3, Figure 4 shows the maximum and minimum principal stresses, σ 1 and σ 3 , at 100% of injection.Figure 4 shows that the maximum of principal stress principally develops zones of tensile stress whereas the minimum of principal stress principally develops zones of compressive stress.For all simulations, the principal zone of tensile stress is located in the host rock directly above the intrusion (Figures 4a,  4c, 4e, 4g, and 4i).The simulation with a bond tensile strength of 15 MPa (i.e., a host rock with a high toughness; Figure 4 e and f) is the only simulation that shows in addition two small tensile zones at the host rock base on each side of the intrusion.In contrast to tensile stress, compressive stress forms two inclined zones, one on each side of the intrusion, which extends laterally away from the intrusion (Figures 4b, 4d, 4f, 4h, and 4j).
The simulation with a bond tensile strength of 1 MPa (i.e., low toughness; Figures 4a and 4b) and the two simulations with an effective modulus of 15 and 30 GPa (i.e., intermediate and high stiffness Figures 4g-4j) display narrow concentrations of high stress where chains of highly stressed particles exist.In contrast, the two simulations with a bond tensile strength of 5 and 15 MPa (i.e., intermediate and high toughness Figures 4c-4f) display a more homogeneous zone of high stresses.
The stresses in the intermediate stiffness (effective modulus = 15 GPa) (Figures 4g and 4h) and higher stiffness (effective modulus = 30 GPa) (Figures 4i and 4j) simulations are very similar in their evolution, patterns, and magnitudes.However, we find that lower stiffness rocks (effective modulus = 3 GPa) (Figures 4c and 4d) are more able to sustain wider and more homogeneous tensile zones compared to host rocks with higher stiffness.In a similar way, but following the opposite trend, the zone of tensile stress becomes wider and more homogeneous when the host rock toughness increases with a bond tensile strength going from 1 to 5 and 15 MPa.For all simulations, the percentage of cracks increases rapidly during the first 30% of injection, regardless of the toughness or stiffness (Figures 5a and 5b).As the injection progresses, we observe sudden increases in the percentage of cracks over short increments of one or two percent of injection.This is particularly evident in the Figure 3.Effect of host rock stiffness on the normalized finite shear strain and spatial crack distribution at 5%, 25%, and 100% of magma injection.The rock stiffness was varied by varying the bond effective modulus, while keeping the bond tensile strength constant at 5 MPa, with a bond cohesion ten times the bond tensile strength.Small black lines perpendicular to bonds represent cracks (broken bonds).Strains on all sub-panels are cut off at 0.05% and normalized.
The total percentage of cracks is lower for a greater rock toughness (Figure 5a).A two-fold increase in bond tensile strength leads to a two-to three-fold increase in the number of cracks.At 100% of injection, a simulation with a bond tensile strength of 10 MPa produces 2.8 times fewer cracks than a simulation with a bond tensile strength of 5 MPa (0.26% and 0.74% of cracks, respectively).The simulation with a bond tensile strength of 20 MPa produces 1.9 times fewer cracks than a simulation with a bond tensile strength of 10 MPa (0.14% and 0.26% of cracks, respectively).Compared to a simulation with an effective modulus of 3 GPa, a simulation with an effective modulus of 15 GPa produces 6.1 times more cracks, whereas a simulation with an effective modulus of 30 GPa only produces 7.0 times more cracks than the 3 GPa simulation (0.74%, 4.48%, and 5.17% of cracks, respectively) (Figure 5b).For all simulations, the lateral extent of vertical surface displacement increases during injection.The sharp increases at 10% of injection correspond to the moment when the subsidence on each side of the uplift surpasses the displacement detection criterion and thus stand out of the noise.We find that host rock toughness does not play a significant role in the surface displacement lateral extent (Figure 5c).The two stiffer host rocks (bond effective modulus of 15 and 30 GPa) show clear stepwise increases in surface displacement lateral extent (Figure 5d).For example, for the stiffer host rock (effective modulus = 30 GPa), there is an increase in the lateral extent of surface displacement from 6.6 to 7.7 km between 23% and 25% of injection.Overall, stiffer host rocks have a narrower extent of vertical surface displacements.
We also observe a stepwise increase in the maximum vertical surface displacement throughout injection for all simulations (Figures 5e and 5f).This can be seen, for example, in the simulations with a bond tensile strength of 20 MPa at 36%-39% of injection, 25 MPa at 48%-53% of injection, and 30 MPa at 59%-62% of injection.Some simulations show sudden and brief decreases up to 0.5 m, for example, between 68% and 72% of injection for a bond tensile strength of 5 MPa and an effective modulus of 3 GPa (Figures 5e  and 5f).
The difference in maximum vertical surface displacement between different values of toughness and stiffness increases throughout injection.At 100% of injection, the difference is about 9 m for toughness variation and 15 m for stiffness variation, which is in both cases less than the largest particle diameter (Figures 5e and 5f).All simulations with a bond tensile strength above 10 MPa have a maximum vertical surface displacement of ≈37 ± 1 m.For the stiffness variation, all simulations with a bond effective modulus ≥10 GPa have a maximum vertical surface displacement of ≈44.0 ± 0.5 m.
In Figure 5, we compared all simulations through the cumulative lateral extent and the maximum of cumulative vertical surface displacement.The temporal evolution of the horizontal and vertical surface displacement profiles can be observed in the movie files provided in Supporting Information S1.

Effect of Source Depth
We simulated several depths of the pre-existing intrusion (0.5, 1.0, and 1.5 km) for two rock types: (a) an intermediate toughness rock, with a bond tensile strength of 5 MPa and an effective modulus of 3 GPa (Figures 6a-6i); and (b) a high toughness rock with a bond tensile strength of 15 MPa and an effective modulus of 3 GPa (Figures 6j-6r).For these simulations, we only illustrate the final state of the injection in Figure 6, that is, when the initial intrusion volume has doubled.
The minimum principal stress in the intermediate toughness host rock shows a slight tensile component with a magnitude below 4 MPa for the three tested depths (Figures 6d-6f).When the source is at 0.5 km depth, the minimum principal stress shows a compressive arched zone with values around 70-80 MPa (Figure 6d).When the source is at 1.5 km depth, two inclined patches of compressive minimum principal stress develop with values above 100 MPa (Figure 6f).The high toughness host rock shows two small zones of tensile minimum principal stress, located within a central 2 km-wide zone for all tested depths (Figures 6m-6o).
After 100% of injection, for a source at a depth of 0.5 km in the intermediate toughness host rock, the host rock is highly fractured (Figure 6g).Multiple tensile fractures open at the surface, and fractures propagate upwards and outwards from the intrusion edges and reach the surface.For a source at 1.0 km depth (Figure 6h), the host rock is still highly fractured and multiple tensile fractures open at the surface.The fractures that propagate upward from the intrusion edges do not reach the surface by 100% of injection.With a source depth of 1.5 km, the normalized finite shear strain is higher than for the shallower sources (Figures 6g and 6i), and the host rock is less extensively fractured (Figure 6i).A single vertical fracture opens at the surface above the center of the intrusion and only small fractures form within its vicinity.
For a source at a depth of 0.5 km in the high toughness host rock, there is very limited fracturing of the host rock (Figure 6p).A vertical fracture opens at the surface above the intrusion center and two fractures propagate upward and outward from the intrusion edges.These two fractures are arched and do not reach the surface.For a source at a depth of 1.0 km, the fracture pattern still consists of a vertical fracture that opens at the surface and two inclined fractures that propagate from the intrusion edges (Figure 6q).However, the two fractures at the intrusion edges are nearly straight.When the source is at a depth of 1.5 km, the normalized finite shear strain is higher than for shallower depths (Figures 6p-6r), but the fracture pattern is similar to the simulation with a 1.0 km source depth (Figure 6r).For all source depths, lower toughness results in a greater degree of fracturing (Figure 7a).The percentage of cracks only shows a dependence to depth for the low toughness rock (bond tensile strength of 5 MPa).For any depth or toughness, the surface displacement lateral extent increases non-linearly through injection (Figure 7b).The lateral extent of surface displacement shows no significant difference between host rocks of different toughness for the same intrusion depth.Deeper sources, however, produce wider surface displacements, and smaller vertical surface displacement maxima (Figures 7b and 7c).

Model Limitations and Advantages
We implement a 2D Discrete Element Method (DEM) application using the software Particle Flow Code (PFC).The designed model setup explores the effect of host rock properties and magma source depth on the deformation and dynamic fracturing induced by the renewed injection of magma into a preexisting intrusion of viscous magma.In DEM models, cracking can occur and coalesce, allowing macro-scale fractures to propagate and affect the future deformation of the host rock.To investigate the effects of the mechanical properties of the host rock, the initial source shape and dimension was constant in all simulations.Furthermore, in all simulations, the injection into the source was performed by applying a constant vertically-upwards rate of 1 m•s 1 to the magma particle inside the conduit.Our model setup represents magma as an incompressible fluid through spherical particles with a high effective modulus but zero friction coefficient or tensile strength.Such assemblage of DEM particles does not precisely represent the flow of liquid magma inside the intrusion from a fluid dynamics point of view, but this approach has been shown to realistically represent the emplacement dynamics of viscous magma domes (e.g., Harnett et al., 2018;Walter et al., 2019).Future work is now required to couple more realistic methods of fluid dynamics modeling with the DEM to simulate the injection of lower-viscosity magma.
We performed one simulation for each of the 13 sets of variables tested (Table 2).An asymmetrical fracture pattern is produced, for example, in the simulation with a bond tensile strength of 5 MPa and an effective modulus of 3 GPa (Figure 3c).Such asymmetry occurs due to the stochastic development of crack distribution in DEM models and may differ from simulation to simulation because of small random differences in the initial particle packing.Three-dimensional FEM models have shown that non axi-symmetrical magmatic pressure sources affect plastic strain patterns in the overburden differently compared to axi-symmetrical sources (e.g., oblate spheroids vs. spheres) (Gerbault et al., 2018).Therefore, we do not expect that the modeled strain, stress, and fracture patterns would be significantly different in 3D models for axi-symmetric reservoirs.3D DEM modeling is now required to investigate the effect of asymmetric reservoirs, but their computation time (days-weeks/simulation) is currently impractical compared to our 2D setup (≈36 hr/simulation).
Our main goal here is to explore the parameter space by varying the (numerical) host rock mechanical properties, rather than to reproduce a particular natural system.We systematically and independently varied the host rock toughness (i.e., the bond tensile strength between 1 and 30 MPa) and the host rock stiffness (i.e., the bond effective modulus between 3 and 30 GPa).Such values correspond to UCS values of 7-203 MPa, tensile strengths of 0.8-21.8MPa and Young's moduli of 2.5-25.3GPa.Such values are within the range of mechanical properties typically measured for intact volcanic rock samples (e.g., Gunzburger & Cornet, 2007;Heap et al., 2020Heap et al., , 2021a;;Heap & Violay, 2021;Perras & Diederichs, 2014;Rajesh Kumar et al., 2011;Yaşar et al., 2010).Those intact properties are typically not representative of largerscale natural systems that contain macro-scale fractures and other discontinuities, and recent work on the cohesion of caldera walls performed with PFC shows that bulk-scale cohesion is likely only a few percent of intact sample strengths (Harnett et al., 2023;Heap et al., 2020).For bond tensile strengths above 5 MPa and bond cohesion above 50 MPa, a single vertical fracture opens at the surface above the intrusion with an aperture on the order of tens of meters.Such a wide surface fracture is not often observed in nature.This corroborates that there may be lower-than-expected bulk strengths, unable to support the height of vertical open surface fractures found in our high toughness simulations.We suggest that the insights obtained on the weaker end-member of our simulated rock strengths are relevant for natural rock masses.
We observed no significant differences in the crack patterns and obtained strength properties in numerical strength tests in which we decreased the mean particle radius from 6.65 to 3.33 and 0.67 m (Figure S2 in Supporting Information S1).On a numerical level, the length of the contact between two particles is controlled by the radius of the two considered particles.The particle radius thus controls the length of a single crack.The particle size does not control how many cracks align to form a fracture.In rock tests, a sample with 6.65 ± 1.5 m particles produces similar fractures as a sample with the same size but with 3.33 ± 0.8 m particles (See Figure S2 in Supporting Information S1).The crack aperture at the particle scale, or the fracture opening at the bulk scale, also (c) as a function of the injection for intermediate toughness (continuous line) and high toughness host rock (dashed line) at three different source depths.The rock toughness did not have a significant impact on the maximum vertical surface displacement or its lateral extent.The depth has a higher impact on these features.The shallower the source, more elevated and narrower the surface displacement.A significant increase in the vertical surface displacement lateral extent as a function of the injection is only visible at shallow depth (black curves).Sudden increases in the percentage of cracks are highlighted by vertical bars, with colors that match the corresponding curves.ten: bond tensile strength; D: magma intrusion depth.
depends on the strength of packing and thus contact properties.This is observable in our results, as all simulations have the same range in particle radius and the produced fractures do not have the same aperture from one simulation to another (Figures S3 to S15 in Supporting Information S1).For simulations with smaller particle radii than used in our simulations, we expect merely a more detailed pattern with smaller-scale fractures in between the large-scale fractures already present.We are confident that the general patterns of stress, strain, fracturing, and response to mechanical property changes highlighted in our work, are independent of particle size.

The Effects of Host Rock Strength on Fracturing
In DEM models, the toughness (compressive and tensile strength in real rock) and the stiffness (Young's modulus) can be independently varied.This helps understanding which physical rock property exerts the most control on damage accumulation and fracture development.For both low toughness (Figure 2c) and high stiffness (Figure 3i) simulations, the host rock is highly fractured following magma injection.However, the low toughness host rock develops wider zones of higher normalized finite shear strain in poorly cracked blocks around which fractures have not yet reached the surface.When the stiffness is varied for a constant toughness, the host rock breaks at different amounts of accumulated strain.For high stiffness host rocks, cracks form at lower values of strain.The normalized finite shear strain in the host rock decreases when new fractures open (compare Figures 2h and 2i).
Both analytical models (with linear elastic behavior) and scaled laboratory experiments (using brittle-elastic gelatine as a crustal analog) have shown that the reduction in lithostatic load as depth decreases affects stress paths and produces linear fractures that rotate and become curved when they approach the free surface (Gudmundsson, 2006;Kervyn et al., 2009;Rivalta et al., 2015).In our simulations, as fractures propagate from the intrusion tips toward the surface, we observe rotation in their orientation which we attribute to the upwards decreasing lithostatic load (Figure 6p).This effect is no longer visible below a depth of 1 km, indicating that the free surface only affects fractures propagating in our simulation in the uppermost kilometer of the crust.
Due to the correlation between Young's modulus and UCS, weak natural volcanic rocks are typically characterized by both a low UCS and a low Young's modulus, whereas strong rocks are characterized by a high UCS and high Young's modulus (see Figure 13b in Heap and Violay (2021)).Figure 8 summarizes that the crack distribution in the host rock varies between two end-members: (a) a highly fractured host rock, with multiple tensile vertical fractures that open at the surface and two well-developed inward-dipping shear fracture zones that develop from the edges of the intrusion; and (b) a poorly fractured host rock, with one vertical central fracture that opens at the surface and propagates downward and, with two inward-dipping fractures that develop from the intrusion edges (Figure 8).The easier the rock is to break, that is, the higher stiffness or the lower toughness, the more the system tends to the highly fractured end-member.
Lower stiffness leads to a less fractured host rock (Figure 3).These two effects oppose each other.There are instances, however, in which volcanic rocks with a low UCS can have a high Young's modulus, such as vesicular (i.e., porous) lavas characterized by a strong and stiff groundmass (Heap & Violay, 2021).If the mechanical properties of the host rocks are known, and the depth and width of the intrusion's roof are known from geophysical data, our model application could provide estimates as to the extent, geometry, and degree of the fracturing within the host rock above or surrounding the intrusion.Additional processes that influence the host rock strength before and during laccolith inflation, such as heterogeneous hydrothermal alteration (Heap et al., 2021b), could furthermore modify the style and amount of fracturing during inflation.
Our simulations allow a better understanding of the dynamics of host rock deformation and uplift patterns observed at geological outcrops as well.At the rhyolitic Sandfell laccolith in East Iceland, the transition from subhorizontal to domed basalt sequences is abrupt and bound by a reverse fault (Mattsson et al., 2018).The radius of the uplifted host rock at Sandfell is 900 m compared to the 700 m radius of the laccolith, and our simulations predict that such inward-dipping reverse faults or shear zones must have grown from the laccolith's edge toward the surface.The uplifted overburden at Sandfell is furthermore disrupted by several normal faults that result in coherent blocks, which would suggest a low toughness host rock and an overall failure pattern more similar to our second, highly-fractured end-member.This is in contrast to the general assumption that basaltic lava flow sequences deform in a stiff, elastic manner above inflating magmatic pressure sources (e.g., Gudmundsson, 2006;Segall, 2010).Such a role of rock stiffness and toughness in controlling fracture patterns was not previously recognized.Galland et al. (2018), for example, only present one conceptual model for the architecture of overburden failure above laccoliths.New field observations should now focus on characterizing host rock strengths and measuring structural indicators of deformation to compare with our model results.
The magma source depth also plays a role in crack distribution, albeit to a lesser degree than the host rock mechanical properties (Figure 6).For host rocks with a low toughness, a deeper source results in less fracturing (Figure 7a) and concentration of cracks in a more localized vertical fracture zone (Figures 6g-6i).This can explain the increase in rock strength with increased lithostatic stress (e.g., Paterson & Wong, 2005) and the resulting depth dependence of contact forces due to gravity in the model.Thus, for a deeper source, the model response changes from the highly fractured endmember to the relatively poorly fractured end-member.Source depth does not affect the crack distribution in a host rock with a high toughness (Figures 6p-6r) or the percentage of cracks (Figure 7a), because the medium is already in the poorly-fractured end-member.
A vertical fracture at the roof center of an intrusion or magma reservoir is generally predicted by analytical models for spherical sources, whereas fractures at the edges are expected for more elongated sources (Browning et al., 2021;Currenti & Williams, 2014;Sigmundsson et al., 2020).This Vshaped fracture zone has been identified in earlier modeling studies of collapse above depleting magma reservoirs (Geyer & Martí, 2014) or pressurization of laccolith-shaped intrusions (Manga & Michaut, 2017).Our model setup produces fractures at the edges and at the top of the intrusion and, even if the source has an elongated shape, fractures develop first at the roof center showing that the fracture initiation is more complex than predicted by linear models.In linear-elastic models, fractures are predicted to initiate at the surface in models that consider the overburden as an elastic plate (Galland & Scheibert, 2013), a very elongated source (Gregg et al., 2012), or by taking into account thermo-mechanical effects within a linear-elastic model (Browning et al., 2021).We observe the development of a near-vertical, through-going fracture that forms by merging a downward-opening tensile surface fracture with a shear band of cracks that propagates upwards from the intrusion's roof.Similar fracture dynamics have been observed in scaled laboratory experiments (Warsitzka et al., 2022).Our results call for a reconsideration of the importance of dynamic fracturing and plastic deformation in numerical models of magma-induced deformation (see also Gerbault et al., 2022).

Implications for Volcano Geodesy and Geo-Resource Exploitation
Our simulations show that the source depth and the host rock stiffness exert a larger influence than host rock toughness on the observed maximum vertical surface displacement (Figures 5 and 7).Host rock toughness does not affect the lateral extent of surface displacement.Host rock stiffness affects the evolution of the maximum surface displacement as injection progresses.Such differences can arise from displacement dissipation in the rock.At shallow depth, there is less attenuation of the displacement in the rock.This attenuation partially arises from the ability of particles to slightly overlap to reproduce elastic deformation (Cundall & Strack, 1979).This is controlled by the bond effective modulus; thus, for stiffer rocks, the amount of displacement that can dissipate in the intact host rock is lower.The difference observed in surface displacement can also arise from the brittle accommodation of deformation in the rock (i.e., fracturing).For any stiffness or high toughness rocks, the rock breaks into large continuous blocks.For example, in Figure 3i, we can identify six coherent, poorly-fractured blocks inside the trapezoidal area defined by the inclined highly fractured zone.In contrast, host rock with a low toughness cannot sustain similar amounts of strain before breaking, and is heavily fractured (Figure 2c).Once the rock is broken into blocks, these blocks are displaced rather than internally strained, resulting in greater surface displacement.By contrast, when the host rock is highly fractured, the stress is accommodated by local shear-mode displacement near the intrusion, resulting in lower surface displacements.In all simulations presented here, we observe abrupt increases in maximum surface displacement and the surface displacement lateral extent over short increments of two or three percent of injection volume.We can relate most of these sudden increments to a sudden increase in the percentage of cracks (Figures 5 and 7).Our model application shows that, even at a constant magma injection rate, the magnitude and/or lateral extent of surface displacement can suddenly change in response to sudden cracking in the host rock.Closer inspection of the strainfracture plots shows that those instances are produced by the sudden coalescence of cracks into macro-scale fractures that cause catastrophic failure of the overburden.In these cases, we observe a transition from cracking near the magma source to direct accommodation of the strain by surface displacement with less cracking within the host rock.This transition from plastic accommodation to elastic bending has also been documented for forced folds above laccoliths and saucer-shaped sills observed in seismic reflection data (Jackson et al., 2013).Further modeling work is now required to investigate the factors that control this transition, including the role of source depth and lithostatic loading.
Our results highlight the importance of considering fracturing to realistically model observed uplift dynamics during sill and laccolith inflation in natural volcanic systems.The non-linearity of cracking observed in our simulations suggests that sudden increases in surface uplift observed by continuous Global Navigation Satellite System (GNSS) measurements in magma-induced deformation events associated with sudden increases in shallow seismic events should be taken as a potential precursor to overburden failure and volcanic eruption.It should be noted that, in our simulations, those abrupt increases in cracking and surface uplift occur without a change in magma injection rate into the inflation source.
Inversions of magma-induced ground deformation using analytical models are sensitive to the estimated magmatic source depth and the simplification of heterogeneous crustal rocks as a homogeneous, isotropic and linear-elastic medium (Masterlark, 2007).Our model application shows that fracture opening affects the relationship between the source depth and the lateral extent and magnitude of the surface displacement.The most significant variations in the lateral extent of the surface displacement occur with source depth variation (Figure 7), in agreement with linear models of deformation (Taylor et al., 2021).However, our results show that the lateral extent of surface displacement slowly increases as injection increases (Figures 5b and 5d).Comparisons between our simulation results and those of linear-elastic analytical and FEM models, will be required to assess the reliability of magma source depth estimates produced using common geodetic models.Our parametric analysis furthermore shows that host rock strength strongly affects the magma-induced deformation patterns and that careful calibration of numerical host rocks to real-world crustal strengths is essential for future modeling efforts.
Our model application shows that an understanding of the mechanical properties of the host rock coupled with an estimation of the source depth can help determine the principal zone of fracturing above a magmatic intrusion.Such fractures can be the preferred path for hydrothermal fluids and ore deposits (Edmonds et al., 2019).Geophysical imaging of the architecture and permeability structure of hydrothermal systems above viscous magma intrusions remains challenging (e.g., Hutchison et al., 2023;Montanari et al., 2017).It is nevertheless essential for the efficient exploration and exploitation of geothermal fields and the interpretation of volcanic unrest events, as exemplified by the hydrothermal unrest that preceded the 2021-2022 eruptive crisis on the Reykjanes peninsula in Iceland (Flóvenz et al., 2022).Our simulations help to locate zones of fracture concentrations above viscous magma intrusions, which form the pathways for hydrothermal fluid migration.Our low toughness, high stiffness end-member, or shallower magma intrusion depths, are more likely to host viable geothermal resources, due to the development of permeability-enhancing pervasive fracturing (Huenges & Ledru, 2011).The localized fracturing in the opposite end-member would favor the concentration of high-grade mineral resources (Rowland & Simmons, 2012).Therefore, if the properties of the host rock above the intrusion and the intrusion source depth are known, our results could be used, alongside other geophysical data, to indicate the best location for the drilling of geothermal wells and the location of high-grade mineral resources.

Conclusions
We present a Discrete Element Method (DEM) model that simulates the mechanical response of the host rock to the inflation of laccoliths by magma injection in the top 1.5 km of the Earth's crust.By systematically varying the host rock stiffness, toughness, and the magma source depth, our results show that the spatial fracture distribution varies between two end-members: (a) a shallow source, or a low toughness or high stiffness host rock that is extensively fractured, with multiple tensile fractures opened at the surface, and inclined fractured zones propagating from the intrusion edges; and (b) a deep source, or a high toughness or low stiffness host rock, that remains relatively unfractured with one central fracture opened at the surface and small fractures at the intrusion edges.The first end-member is favored in the case of a high rock stiffness, a low rock toughness, or a shallower intrusion.The estimation of host rock mechanical properties and source depth can thus help to determine the preferred path for hydrothermal fluids or ore deposits.Abrupt increases in the maximum and lateral extent of vertical surface displacement can arise due to the fracturing of the host rock, even at constant magma injection rates.As a result, rapid increases in surface displacements measured by geodetic monitoring equipment above laccolith intrusions may not be related to increases in magma injection rates and, if contemporaneous with local seismicity, may rather be the result of host rock fracturing that could lead to catastrophic overburden failure and volcanic eruption.Our model application can act as an important tool for assessing a more realistic response of a heterogeneous crust to magmatic intrusion, as well as the architecture of hydrothermal systems or mineral deposits above shallow magma intrusions.

Figure 1 .
Figure 1.General geometry of the 2D model setup before injection.(a) The layers are for visualization purposes only and do not correspond to different mechanical properties.Each layer corresponds to 200 m except the lower one which is 100 m.(b) Zoom on the laccolith at the bottom of the model.The blue part of the box corresponds to a permeable wall that only allows particles to go through in the upward direction.

y
refers to the maximum vertical surface displacement and U y to the surface displacement.a Mean value of 15 rock tests.b 1st end-member: Highly fractured pattern, similar to Figure 2c or Figure 3i.2nd end-member: Poorly fractured pattern, similar to Figure 2i.Mixed: Fracture pattern not clearly in an end member pattern, similar to Figure 2f. to those outliers in all shown stresses.Stresses were cut off at ≈0.5% of positive and negative values, meaning 15 MPa for positive values and 150 MPa for negative values.

Figure 2 .
Figure2.Effect of host rock toughness on the normalized finite shear strain and fracture pattern at 5%, 25%, and 100% of magma injection.The rock toughness was varied by varying the bond tensile strength, while the bond effective modulus was kept constant at 3 GPa for all three simulations shown here.Small black lines perpendicular to broken bonds represent cracks (broken bonds).Strains on all sub-panels are cut off at 0.05% and normalized.

Figure 4 .
Figure 4. Effect of the variation of host rock toughness and stiffness on the maximum and minimum principal stresses σ 1 and σ 3 .Compression is negative and tension is positive.The stress amplitudes are all cut off at the same value to avoid a flattening of the stress patterns and allow comparison between simulations.In all simulations the bond cohesion is ten times larger than the bond tensile strength.

Figure 5 .
Figure 5. Evolution of the percentage of cracks (a-b), the lateral extent of the vertical surface displacement U y (c-d), and the maximum vertical surface displacement U max y (e-f) as a function of the amount of injection for host rock toughness variation (a-c-e) and stiffness variation (b-d-f).Sudden increases in the percentage of cracks are highlighted by vertical bars, with colors matching the corresponding curves.ten: bond tensile strength; E*: effective modulus.

Figure 6 .
Figure 6.Maximum and minimum shear stress σ 1 , σ 3 and normalized finite shear strain at 100% of injection as a function of the source depth for intermediate toughness (a to i) and high toughness (j to r) rocks.(a to f and j to o) Orange-colored patches indicate compression and blue-colored patches indicate tension, color bars are cut off at the same value.(g to i and p to r) Small black lines perpendicular to broken bonds represent cracks.All strains are cut off at 0.05 and normalized.In all represented simulations the effective modulus is 3 GPa.

Figure 7 .
Figure 7. Evolution of the percentage of cracks (a), the lateral extent of the vertical surface displacement (b) and the maximum vertical surface displacement U max y

Figure 8 .
Figure 8. Representation of the two end-members between which the fracture pattern varies as the toughness, stiffness or source depth vary.Note that horizontal layering in the particle packing is only for visual purposes and does not represent changes in mechanical properties.

Table 2
Summary of the Key Results for the Thirteen Simulations Where the Source Depth, the Effective Modulus E*, the Bond Tensile Strength ten, and the Bond Cohesion coh Were Varied