Mechanical Modeling of Pre-Eruptive Magma Propagation Scenarios at Calderas

. The simplest two-dimensional (2D) trajectory models are streamlines perpendicular to the least-compressive stress Abstract Simulating magma propagation pathways requires both a well-calibrated model for the stress state of the volcano and models for dike advance within such a stress field. Here, we establish a framework for calculating computationally efficient and flexible magma propagation scenarios in the presence of caldera structures. We first develop a three-dimensional (3D) numerical model for the stress state at volcanoes with mild topography, including the stress induced by surface loads and unloading due to the formation of caldera depressions. Then, we introduce a new, simplified 3D model of dike propagation. Such a model captures the complexity of 3D magma trajectories with low running time, and can backtrack dikes from a vent to the magma storage region. We compare the new dike propagation model to a previously published 3D model. Finally, we employ the simplified model to produce shallow dike propagation scenarios for a set of synthetic caldera settings with increasingly complex topographies. The resulting synthetic magma pathways and eruptive vent locations broadly reproduce the variability observed in natural calderas.

Extending the stress calibration strategy by Rivalta et al. (2019) to 3D would pave the way to forecast dike pathways in 3D at any arbitrary volcano. A preliminary step is to set up 3D stress and dike trajectory models that are computationally efficient for the large number of simulations needed by the stress calibration procedure. In this study, we first develop computationally efficient 3D stress field calculations for scenarios with topographic reliefs. Then, we develop a fast, semi-analytical 3D dike propagation model that approximates the sophisticated model by Davis et al. (2020Davis et al. ( , 2021 but retains the simplicity of 2D streamlines and can also backtrack a dike trajectory from eruptive vent to magma chamber. Finally, we show how to integrate all these models to produce realistic pre-eruptive magma propagation scenarios. We focus on calderas, setting up synthetic topographies inspired by natural systems.

Method Formulation
We assume a homogeneous, isotropic and linearly elastic medium as the host rock, described by rock density ρ r , Young's modulus E and Poisson's ratio ν. g is the acceleration due to gravity. Symbols and parameters are defined in Table 1.

A Modular Approach to Understanding Stress States
We describe the state of stress within the host rock by a stress tensor σ ij . Tensional stresses are positive. σ ij is diagonalized to retrieve magnitudes, σ 1 , σ 2 , σ 3 , from most compressive to least compressive, respectively, and eigenvectors, 1 , 2 , 3 , which identify the orientations of the principal stress axes.
We build our 3D stress model following the first-order linear approach by Rivalta et al. (2019), who expressed the elastic stress field σ ij of a volcanic region as the superposition of perturbations from a background stress state 0 , each stemming from a different stress-generating mechanism. The approach neglects coupling between the stress sources. We limit our analysis to tectonic stresses and gravitational loading/unloading because dike patterns can often be explained by a combination of the two mechanisms (Corbi et al., 2015;Heimisson et al., 2015;Maccaferri et al., 2017;Neri et al., 2018;Roman & Jaupart, 2014). Not including other mechanisms, such as pressurized magma reservoirs or faults, has the advantage of limiting the number of parameters in the model, while retaining the stress mechanisms with the largest influence. More contributions can be easily added, if needed in specific cases.
We write the stress tensor at any point in the crust as: where the terms on the right side arise, respectively, from the regional tectonic stress (T) and the gravitational loading/unloading (G).

10.1029/2022JB025956
3 of 19 The first step is to define the unperturbed state of stress, 0 , before any of the sources on the right hand side of Equation 1 became active. There are two main assumptions in literature: a laterally-confined medium, that is, no lateral strain can be produced after gravity is turned on (e.g., Martel & Muller, 2000;Savage et al., 1985), resulting in a vertical 1 : or a lithostatic stress state: Field measurements of subsurface stress (Jaeger et al., 2007) lie somewhat in between those two assumptions. Therefore, 0 can be written as: where ∈ [ (1− ) , 1 ] (Jaeger et al., 2007;J. R. Muller et al., 2001;Slim et al., 2015). In this study, we set k = 1 and assume a lithostatic unperturbed stress.
The second step is to superimpose the tectonic stress, expressed in terms of three independent components , , , here assumed uniform (e.g., McKenzie, 1978;Müller et al., 1992).
The third step is to consider gravitational stresses associated to surface loading or unloading. This has often been modeled by distributions of normal forces onto a half-space (Dahm, 2000b;Maccaferri et al., 2014;Neri et al., 2018), which, however, neglect the shear stresses imposed by the topography and provide no information on the stress within the topography itself (McTigue & Mei, 1981). More sophisticated analytical solutions exist,  & Mei, 1981;Savage et al., 1985) or only for simple topographies (McTigue & Mei, 1987). Stress due to surface loading/unloading decays over a vertical distance that scales with the radius of the topographic feature (e.g., Jaeger et al., 2007;Pollard et al., 2005;Roman & Jaupart, 2014). Consequently, principal stresses can change in both intensity and orientation over short distances. This has several implications discussed later (Section 2.2.1). Martel and Muller (2000) and Slim et al. (2015) described how to implement topographic loads within Boundary Element (BE) models, where the topography is discretized into a mesh of dislocations. They considered the effect of topographic loading as akin to cutting an infinite body subject to gravity in two halves along a surface defined by the topography. The gravitational stress imposed by the upper half onto the lower one is then subtracted from the background stress of the body (Martel & Muller, 2000, Figure 3). In practice, this is achieved through imposing boundary conditions on the BEs, depending on the coordinate z of their midpoints and the rock density, which control the overburden or excavation pressure imposed by the topography.
One important point in models such as Martel and Muller (2000) is that the boundary conditions at the BEs representing the topography are univocally fixed only once the datum level, that is the unperturbed surface before any topography is created, is set. This was rarely clarified in past applications (e.g., Chadwick & Dieterich, 1995;Neri et al., 2018;Urbani et al., 2017). Identifying such surface is not always trivial but critical, as different choices lead to different outcomes for the displacement and stress field. We show this in Figure 1a, where we compare 1 from the analytical solution by McTigue and Mei (1981) for a valley adjacent to a ridge under plane strain assumption to 2D numerical models where the datum level is set to, successively, the flat extremes of the profile, the ridge summit and the valley bottom. The first model shares the same assumption on the datum level with the analytical solution, hence the good agreement for that case. Such assumption is straightforward to adopt when the topography becomes uniformly flat away from the loaded/unloaded region. However, this is not always the case, and the optimal choice of datum level may depend on the situation. Take, for example, a caldera lying on a coastline, which divides two regions, the mainland and the sea floor, at different elevations. We consider a similar case in our synthetic scenarios, and we solve the ambiguity in the datum level by setting it to the ground elevation before the caldera was formed: this coincides with the sea level in that case. If, for instance, we were to study the formation of an edifice and, later, of a caldera at its summit, we would first set the edifice datum level at its base, and then set the caldera datum level at the edifice summit. Consequently, the topography preceding the reference event (in our scenarios, the caldera formation) informs the datum level.
A further issue regarding the calculation of surface loading/unloading stresses is that they are not immutable. Volcanic regions host a variety of stress-generating and stress-relieving mechanisms acting on different time scales. For example, the build-up of a volcanic edifice consists of progressive accumulation of eruptive material that loads and stresses the underlying crust (McGuire & Pullen, 1989), while, at the same time, magmatic intrusions, earthquakes and inelastic processes tend to relax shear stresses and homogenize principal stresses (Chadwick & Dieterich, 1995). Quantifying stresses within large topographic loads at a particular point in time is thus non-trivial. Here we avoid this issue by focusing on calderas that we assume have formed relatively recently in the history of the volcano, and consider otherwise only mild topographies, so that modeling dike propagation within edifices is not necessary. We elaborate further on this point in Section 4. We note that when we use the term "caldera," we are referring to the general surface depression that is associated with all calderas. Differences in the origin, structure and setting of calderas (e.g., Acocella, 2007;Cole et al., 2005) are neglected.
We compute ( ) in Equation 1 following Martel and Muller (2000), Slim et al. (2015). We employ the 3D BE tool Cut & Displace (Davis et al., 2017, based on the displacement discontinuity method by Crouch et al. (1983). We use DistMesh by Persson and Strang (2004) to discretize the topography into a mesh of triangular dislocations (Nikkhoo & Walter, 2015), acting as BEs. The 3D mesh needs to be larger than the region of interest, so that its edges are distant enough from the volume where we compute the stress. We find that a mesh with a diameter three times the lateral extent of the studied region is enough for that purpose, and we adopt this choice in all our models. If a coastline is present, the outer mesh tapers to two horizontal surfaces at different height, representing the far-field mainland and the far-field sea floor. Once the datum level is fixed, stress boundary conditions are imposed on each BE as previously described. The load imposed by the water column on the bathymetry is also included.
Calderas are usually filled with eruptive material or sediments over time (e.g., Hildreth et al., 2017;Orsi et al., 1996). Our model can account for this in several ways: the original buried caldera floor may be meshed as the reference topographic relief, and the corresponding BEs may be loaded accounting for the pressure deficit 6 of 19 due to the density contrast between the deeper host rock and the layers above. Alternatively, the current caldera topography may be meshed as the reference topographic relief, and the unloading pressure resulting from the missing mass due to lower density infill is factored in the boundary conditions. Calculations for these options for a synthetic caldera (Figure 1c) show good agreement except in the proximity of the caldera rim. Here we follow the former approach in one scenario, as illustrated later.
We remark that some of the stress sources we neglect, such as magma reservoirs, are in principle straightforward to include in our BE model. In order to show the minor relative influence of such sources, we compare in Figure 1b the orientation of 3 for three different models: one without and two with a pressurized, spherical magma chamber, with overpressure of 10 and 100 MPa, all involving the same surface unloading and tectonic stress. Only with extremely large overpressures the effects of the pressurization are felt at a distance of up to one source diameter. This validates in 3D a similar argument by Rivalta et al. (2019) (see their Figure 1).

Simplified Analytical Model (SAM)
Next, we develop a computationally-efficient 3D dike propagation model that provides a 3D equivalent to 2D 3 -perpendicular streamlines. There is no straightforward method to compute streamlines in 3D, as the direction of 3 alone identifies a surface, while the direction of propagation on that surface remains undetermined. Davis et al. (2020Davis et al. ( , 2021) developed a point-wise, analytical dike trajectory calculator, similar to Sigmundsson et al. (2015) but fully 3D and more comprehensive in terms of factors considered. Its purpose was to justify why an observed dike took a specific direction depending on the magma buoyancy and the external state of stress, and falls short of being a propagation model. Here we turn that approach into a simplified 3D propagation model that can also backtrack dike trajectories downward from a vent to the magma storage region. We henceforth refer to our model as the "Simplified Analytical Model" (SAM).
In the analytical model by Davis et al. (2020Davis et al. ( , 2021, propagation of the tip-line of a dike occurs when the local mode I stress intensity factor, K, is larger than the fracture toughness, K C , of the host rock (e.g., Secor & Pollard, 1975). The dike is represented as a tensile penny-shaped crack with a fixed volume, V, and radius, c. It is assumed that external stress varies linearly in every direction over the crack surface, and that internal pressure varies linearly with z proportional to ρ m g sin β, where β is the crack dip. In such case, K can be written as: (Tada et al., 2000), where Δγ max is the maximum value over all orientations across the crack plane of the "pressure gradient", Δγ, calculated as the difference of the external stress and internal magma pressure over the crack diameter, and α is the angle spanning the circumference of the crack away from the direction of Δγ max (see Figure 2b). The second contribution in Equation 5, which is largest for α = 0, determines the maximum of K and, thus, the direction of propagation of the crack. If R K = K/K C > 1, the crack propagates (see Figure 1 in Davis et al., 2020).
In SAM, we simplify such approach by forcing the dike to open against the local 3 , and calculating K simply as This is equivalent to neglecting the role played by the dike volume and K C in determining whether the dike will advance. On the other hand, the buoyancy force contributes to Δγ, and plays a role in determining the direction of propagation on the 3 -perpendicular surface.
In a Cartesian reference frame, where the z-axis is positive upward (Figure 2a), we calculate forward dike trajectories (FTs) as "paths of local steepest ascent," corresponding to the steepest increase of Δγ, as follows: Evolving topography: a 1-km-deep axisymmetric caldera is refilled by 1/3 of its original depth. 3 orientation and topographic profiles for two mechanically-equivalent models of caldera unloading with different reference topographic relief and boundary conditions. (c) Importance of reservoir: 3 orientation for three models involving a 1-km-deep axisymmetric caldera and vanishing tectonic stress. Two models include a 6-km-deep spherical magma reservoir of 1.5 km radius with overpressures ΔP = 10 MPa (red) and 100 MPa (green) respectively; one has no reservoir (blue).
7 of 19 1. We produce a stress model for the hosting medium (Section 2.1). 2. We choose a starting point F 0 for the dike (for instance, at the edge of a magma reservoir). 3. We compute σ 3 and 3 at F 0 and identify the local surface Σ perpendicular to 3 . The dike is then defined as a penny-shaped crack of radius c lying on Σ (Figure 2a). 4. We generate a ring of n regularly-spaced observation points O i , i = 1…n along the dike tip-line ( Figure 2b). 5. We calculate 3 at each O i and use it to calculate Δγ for every point on the dike tip-line as: where , are the vertical coordinates of points O i , O j , with O j antipodal to O i . 6. We calculate K i at each O i according to Equation 6 and determine the point F 1 where K i = K max . This will identify the direction of propagation of the dike ( Figure 2b). Such direction coincides with that of the maximum pressure gradient across the plane of the crack. Note that negative K are always predicted at some O i and imply unrealistic interpenetration of the crack faces. This poses no issue, however, since we are only interested in finding K max .
We reiterate the previous steps taking F 1 as the current F 0 and produce a chain of points identifying the trajectory of the dike. The dike stops once at least one of the observation points generated in step 3 reaches a minimum distance threshold (MDT) between the observation points and the mesh, in order to prevent artifacts or singularities in the stress calculations. This is a characteristic issue of BE models, and can be mitigated with finer meshing (Slim et al., 2015). Here we fix the MDT to 800 m away from the nearest BE, as this is the average size of the dislocations of the mesh we employ. Dikes may be propagated past their F A until they hit the surface at a "projected" arrival point, , assuming that they maintain the dip and strike calculated at F A ( Figure 2b). This is akin to assuming that dikes do not have the space to adjust to the local stress field in the last ∼1 km before reaching the free surface. Moreover, a SAM dike is forced to stop if the trajectory becomes horizontal, or if the difference in the strike and dip angles between the current direction of propagation and the one at the previous step is larger than a given threshold. This prevents abrupt turning of the dike pathways.
SAM trajectories depend on two parameters, c and n. We found that values of n equal or greater than 12 lead to nearly identical dike pathways; we set n to 12 in all scenarios calculated later. In contrast, different c lead to different trajectories and arrival points for the same starting points and stress field. Large c (e.g., >2 km if the dike starting point is 10 km deep) sample the stress field in too few points and approximate Δγ too coarsely to produce accurate trajectories, while very small c (e.g., <50 m for the starting depth mentioned above) follow principal stress directions nearly point-wise, but are more computationally expensive. We show later how c may be calibrated to better match a more sophisticated dike propagation model. SAM also allows for the propagation of anti-buoyant dikes, that is, dikes filled with ρ m > ρ r propagating downward through the crust. Dike trajectories, however, cannot be backtracked by simply inverting the density contrast between magma and rocks: an anti-buoyant dike starting from the arrival point of a buoyant one and propagating downward with the same c and n will not pass through the same points (see Figure 2c), even if the difference between forward trajectories (FTs) and backtracked trajectories (BTs) decreases for smaller values of c.
We backtrack FTs, from known arrival points F A and with assumed parameters c B and n B , as follows: 1. Starting from F A , we find a candidate point B C at a distance c B such that the scalar product between 3 at B C and the vector pointing from B C to F A is minimal (Figure 2d). 2. We run one step of the forward model from B C and calculate the vector between the predicted and actual F A ; we then shift B C by that same vector and iterate this procedure until the desired precision is attained. B C is taken as the first point B 1 of the BT. 3. The algorithm stops as soon as a specific requirement is satisfied: for instance, the current B j falls within the known magma storage region. The lastly-recovered point of the BT becomes then the "backtracked starting point" (BSP) (Figure 2d).
The first step of the algorithm is modified when starting from a point lying on the free surface, as we no longer fix the distance between B C and to a specific c B , but let it vary over a specific range (for a FT with given c, we find a 0-3c range enough for our purpose).
We tested the method against known FTs, and found that it is able to retrieve each F 0 within a range of a few tens of meters (∼0.2%-0.5% of a 6-km caldera radius) if starting from F A , and a few hundreds (∼2%-5% of the same caldera radius) if starting from , provided the same radius c of the forward model is employed (c B = c). If that is not the case, the distance between actual and BSP (Δ BSP = |F 0 − BSP|) increases with the difference between the backtrack radius c B and c.

Three-Dimensional Intrusion Model (TIM)
We later validate SAM against the full-3D numerical dike propagation model by Davis et al. (2020) and Davis et al. (2021). The model needs the dike volume (V), assumed constant during the propagation. Here, the dike starts as a penny-shaped crack centered at a specific starting point and arranged according to a starting dip and strike; these can be either arbitrary or coincide with the local 3 . The dike starting radius is taken as 0 = √ ∕1.6 . The dike is meshed, and R K is computed at every tip-line BE ; the tip-line is advanced or retreated by an amount proportional to the local R K , depending on its sign, and the crack is remeshed. The crack can also bend out of its plane according to the maximum circumferential stress criterion Pollard et al., 2005). The dike can thus advance along complex trajectories and change its shape in the process. We refer to this model as "Three-dimensional Intrusion Model" (TIM).
TIM relies on finer discretization (at the scale of individual BEs) when calculating K. Comparing the two models is therefore critical to verify the validity of the approximations in SAM, especially at shallow depths, where even minor topographic features have a non-negligible influence and lead to more heterogeneous stress gradients (see Section 2.1).
Before comparing TIM and SAM trajectories, we illustrate how to combine the stress and dike models introduced so far into synthetic scenarios of dike propagation.

Configuration of the Dike Propagation Scenarios
We produce a total of nine synthetic scenarios (Tables 2 and 3). We first generate a stress model, evaluating which stress mechanisms are most relevant. Here, as discussed in Section 2.1, we limit our analysis to tectonic stresses and gravitational loading/unloading.
We consider increasingly complex topographies with a caldera located at the origin of the Cartesian reference frame (see Figure 2a). We employ four main topographic settings, each used in one or more scenarios: • Setting 1: a flat topography with a circular caldera of radius R C = 6 km and maximum depth d = 500 m. The depth of the caldera, which has steep slopes and a flat floor, varies with r according to: • Setting 2: we add a coastline, modeled as a steep elevation step along the y-axis. In this way, we break the axial symmetry of the previous setting. The bathymetry lies 100 m below the datum level. The caldera has R C = 6 km, d = 450 m, and depth varying with r as in (1).
10.1029/2022JB025956 9 of 19 • Setting 3: we maintain the bathymetry of (2), but we include two hills (heights 791 and 355 m, base diameter ∼15 km). The caldera has R C = 6 km and d = 424 m. The caldera shape is made irregular by adding Gaussian noise to Equation 8. In one scenario we model a topography evolving from (3) to (3b), where the caldera is partially refilled, its maximum depth changing to d = 221 m. This setting is inspired by the morphology of Campi Flegrei caldera. • Setting 4: an elliptic caldera with d = 150 m, semi-major and semi-minor axes a C = 8 km and b C = 4 km, respectively. A circular resurgent dome with h = 150 m and 4.8 km diameter is located 3 km offset from the caldera center. The external topography has some gently-sloping hills (the maximum height is 157 m), but no bathymetry. This setting is inspired by the morphology of Long Valley caldera.
We calculate the gravitational loading/unloading as described in Section 2.1, using E = 15 GPa, ν = 0.25 and setting ρ r as in Table 3. Then, we superimpose to the resulting stress field the tectonic stress components , different for each scenario.
Next, we choose a model of dike propagation and define the needed input. TIM needs dike volumes (V k ), magma densities ( ) , K C of the host rock and a starting geometry for the kth dike. We use K C = 70 MPa⋅m 1/2 . Starting dike radii ( 0 ) are determined by V k (see Section 2.2.2 and Davis et al., 2021). SAM needs c and ρ m .
We use the first three scenarios, "Vertical-TIM," "Lateral-Dike," and "Complex-Coastline," to compare the performance and features of TIM and SAM. In the additional six scenarios, we produce only SAM dike pathways with fixed c = 1.2 km and ρ m = 2,300 kg/m 3 . We start with the most simplified topography ("Circular-Caldera"). Then, we progressively add new elements, such as a coastline ("Simplified-Coastline", "Tectonic-Shear"), hills and caldera refilling ("Refilling-Caldera", "Two-Reservoirs") and a resurgent dome ("Elliptic-Caldera"), studying their impact on dike trajectories. In three scenarios, the arrival points of SAM dikes are projected past the MDT to the free surface (see Section 2.2.1). All scenarios involve tensional stresses, whose principal axes coincide with the coordinate axes except for Tectonic-Shear (Table 2).
We fix the location, depth (z res ) and radius (r res ) of the magma reservoirs, which constitute the rock volumes where dikes depart from. We remark  that here the reservoirs have no contribution to the stress field. All magma reservoirs are circular, sill-like and centered at the origin of the reference frame. In Elliptic-Caldera, however, we consider a vertically-elongated reservoir centered below the summit of the resurgent dome.
The number of simulated dikes (N) varies among the scenarios (Table 2). Dike starting points are described by depth 0 , radius 0 = res and angle 0 , k = 1,…,N, according to the cylindrical reference frame in Figure 2a. In the most simplified scenarios, we assume equally-spaced starting points for dikes. In the most complex scenarios, we randomize the starting points by drawing 0 from an uniform distribution. Starting depths coincide with the depth of the magma reservoir ( 0 = res ) , with two exceptions. In Two-Reservoirs, we consider two different starting depths, with the aim of representing two separate magma storage volumes. In Elliptic-Caldera, we draw 0 from a Beta distribution (e.g., Johnson et al., 1994) skewed toward the top of the reservoir (see Figure 5f). This is to simulate a case where dike nucleation probability may change with depth.

SAM and TIM Comparison
We now proceed to validate SAM against TIM to assess under which conditions the two models produce matching dike pathways. We use Vertical-TIM, which offers the simplest topography, and Lateral-Dike, which offers the most complex one, to compare TIM and SAM pathways under different starting conditions and settings. Then, we use Complex-Coastline to calibrate c in SAM.
If TIM dikes start misoriented with respect to the external stress field, they can progressively adjust to it as they advance, while SAM dikes start and remain perpendicular to 3 . This can lead to discrepancies between SAM and TIM dike pathways. We show this in Vertical-TIM (Figure 3a), where three vertically-oriented TIM dikes with different volumes (V k ) and starting radii and two SAM dikes with different c propagate from the same starting point and with the same ρ m (Table 3). In Figure 3a, TIM and SAM dikes first diverge, and later become roughly parallel, as TIM dikes adjust to the stress directions. Dikes with larger volumes require larger distances to do so, as already captured by 2D models (Dahm, 2000a;Maccaferri et al., 2010Maccaferri et al., , 2019. We also notice how the SAM dike with the smallest c follow the stress field more closely. In Lateral-Dike we show a situation where SAM captures 3D propagation as well as TIM. We run a TIM dike starting beneath a topographic high, and compare it to a SAM dike with radius c = 1.2 km starting from the same point. In this model, we set both dikes to be weakly buoyant (ρ r − ρ m = 100 kg/m 3 ) and start aligned to the local stress directions. In these conditions (Tables 2 and 3), they both propagate laterally along similar trajectories, as dictated by the external stress and the low magma buoyancy: such behavior may not be captured by 2D dike models. In Figures 3e and 3f we observe that K values in SAM can be very different from the ones in TIM, and the SAM dike follows a longer, zigzagging pathway. This is due to the large c employed, which makes the dike advance too far to capture at each step the heterogeneity of the pressure gradient. Notwithstanding these differences, the overall directions of the pressure gradient (orange arrows in Figures 3e and 3f) are consistent, and the dikes follow each other closely even at shallow depths. In a test not reported here, we run the same scenario with a larger magma buoyancy (ρ m = 2,300 kg/m 3 ), and both TIM and SAM dikes ascended toward the free surface instead of propagating laterally. This shows how accounting for the magma buoyancy force in SAM makes it different from a simple "3D streamline" approach, as SAM dikes do not necessarily follow 1 .
In Complex-Coastline (Figures 4a-4d), we study a case where TIM dikes start optimally-oriented (i.e., perpendicular to 3 ). We run nine TIM dikes with different V k , 0 and ρ m (Table 3), and compare them to forward SAM trajectories. Despite the V k , 0 and being different from one dike to another, the arrival points and final orientations of the SAM dikes are consistent with the outcomes of the TIM dikes, and SAM trajectories follow closely TIM ones. Such match is closest when we take =̄0 , that is, the average of the 0 (Figure 4c). In order to refine our calibration of c, we perform a comparison between TIM dikes and backtracked SAM dikes, evaluating how accurately their starting points are recovered with different values of backtrack radius c B (Figures 4e and 4f). We find that the performance of our backtracking method in recovering the SP of the TIM dikes depends on the c B we employ (see Table 1 for abbreviations). Both large (>1.2 km) and small (<0.6 km) c B perform poorly. On the other hand, the distance between SP and BSP of each dike, Δ BSP , is smallest for c B equal or close to 0 = 880 m (black vertical line in Figure 4e). The minimum of Δ BSP for all dikes except for the one with the smallest V k ( leading to the most accurate BSP for the kth dike, shows that the best-fit line comes close to the bisector of the quadrant and, thus, = 0 = √ ∕1.6 provides a good estimate for the optimal radius c in SAM (Figure 4f).
In summary, SAM provides trajectories close to TIM dike trajectories only when the latter are well-oriented within the external stress field. In that case, the two models compare well even if the predicted values of K are very different. The optimal c for SAM may be chosen on the basis of the volumes of TIM dikes. The implication  Figure. (c) NW-SE view of (b) looking from the direction shown in (b) as an orange arrow. TIM dike represented as superposition of red meshes from five steps in the dike simulation, from start to end. Each step of SAM pathway is a green circle. (d) Outlines of the five steps of TIM pathway shown in (c). SAM cracks are superposed in gray. (e) Values of K computed along tip-line of TIM meshes, as well as K gradient directions for each step (black) and K gradient direction averaged over whole pathway (orange).
Step 1 in (d) not shown here. (f) Values of K computed at observation points along tip-line of SAM cracks, as well as K gradient directions for each step (black) and K gradient direction averaged over whole pathway except for last step where dike stops (orange). For host rock and magma properties see Table 3.

of 19
is that, in a real scenario, knowledge on the volume of actual dikes could inform the choice of c for both forward and backward SAM. We add that, in Lateral-Dike, the running time of one step of SAM is ∼100 times faster than that of one step of TIM.

Results
In the simplest model (Circular-Caldera, Figure 5a), dike trajectories are deflected by the gravitational unloading associated to the caldera, and their arrival points punctuate its rim. The tectonic extension is higher along the x-axis, and this leads to the spacing between neighboring arrival points becoming smaller when closer to that axis, even if the starting points are equally spaced.
In Simplified-Coastline, the presence of a coastline between two flat regions at different heights has an evident impact on dike trajectories, which are still deflected away from the caldera, but end up mostly on the mainland  Table 1). Black dotted line marks the average of 0 of TIM dikes. Colors are the same of TIM and SAM arrival points in (b) and (c), and numbers in the inset follow the order of Table 3. (f) : c B yielding the smallest Δ BSP versus starting dike radius for each dike. The red line fitting the data is compared to the bisector (blue line).

of 19
( Figure 5b). Only the dike starting farthest away from the mainland manages to reach the sea floor. In particular, there is a concentration of arrival points close to the coastline. The effect of deviatoric tectonic stress is most apparent in Tectonic-Shear (Figure 5c). Here, the least-compressive principal tectonic stress axis roughly strikes along the bisector of the second and fourth quadrants (NW-SE). Arrival points cluster about such axis, both on the mainland and on the sea floor. Caldera refilling and the presence of a resurgent dome cause an inward shift of dike trajectories. In Figure 5d (Refilling-Caldera), green dikes are still deflected by the caldera unloading, but all reach the surface along or within the caldera rim, some ending up on the resurgent dome. Topographic loads outside the caldera tend to attract dikes from both red and green sets.
Dikes departing from deeper storage regions, as in Two-Reservoirs (Figure 5e) show the same pattern as in the previous scenarios, reaching the surface farther away from the caldera.
Dikes in Elliptic-Caldera (Figure 5f) feel the competing influence of the elliptic caldera and the loading due to the resurgent dome and the hill west of the caldera. The synthetic vents cluster in two areas, the larger adjacent to the dome and the minor close to the caldera rim and the hill. No vents are present at the top of the dome.
In most scenarios, many dikes stop before reaching the MDT (Table 2) when the interplay between the buoyancy force and the external stress gradients is no longer sufficient to drive the dike upward. Dike arrest is often associated to gravitational loading (topographic highs): in Refilling-Caldera and Two-Reservoirs, most dikes ascending below the highest hills stop before reaching the MDT. This is consistent with the outcome of Lateral-Dike (Figures 3c and 3d), where both TIM and SAM dikes stop ascending and propagate laterally beneath a topographic load before stopping.
In summary, topography plays a dominant role in controlling dike pathways in our scenarios. Even relatively small topographic features, such as the ∼5-km-wide resurgent dome in Elliptic-Caldera (Figure 5f), influence close trajectories over a distance comparable to their width. In all scenarios, dikes are consistently deflected away from surface unloading and attracted by surface loading. Tectonic stress also influences dike orientation and clustering of arrival points, with a more evident impact in the simplest scenarios (Circular-Caldera, Simplified-Coastline, Tectonic-Shear).

Discussion and Conclusions
We have shown how our newly-developed "elementary" dike propagation model (SAM) well reproduces trajectories calculated with a sophisticated numerical model (TIM) by Davis et al. (2020); Davis et al. (2021) (Figures 3b,3c,3d, and 4), and can effectively model 3D dike pathways in synthetic calderas with tectonic stress and mild surface loading/unloading ( Figure 5). In particular, SAM and TIM trajectories are similar if TIM dikes start optimally-oriented to the external principal stress directions (Figure 3a), since SAM dikes are always oriented perpendicularly to the local 3 . Moreover, if stresses change over a distance smaller than c, the calculation of the pressure gradient (Section 2.2.1) and, consequently, SAM trajectories, will be more approximated. Large c values, however, are still reliable in our scenarios, since loads/unloads with large horizontal extent (>5 km) cause smoothly-changing stresses at the scale of most SAM dikes shown here (see Section 2.1). Loads/unloads of small extent (<1 km) would cause rapidly-changing stresses at that same scale, but their effect is significant only at shallow depths and can be neglected here, as we stop dikes at the MDT (Section 2.2.1). In this regard, fixing the MDT determines what topographic details are worth considering in our models. Dike propagation in both models is controlled not only by the gradients of external stress, but also by magma buoyancy. SAM is also able to backtrack dike trajectories from a vent to the magma storage region.
Due to our simplifying assumptions, our models have many potential limitations. The assumptions include homogeneous elastic parameters for the host rock. Rigidity and density layering may substantially affect dike propagation. For instance, dike trajectories can be deflected when crossing interfaces between layers with strong rigidity contrasts, as shown in 2D by Maccaferri et al. (2010). Further studies are needed to grasp the effects of layer interfaces in 3D. Nevertheless, as shown by Mantiloni et al. (2021) through analog experiments, homogeneous models well reproduce the observed pathways provided that "effective" stress parameters are employed, rather than those actually imposed on the gelatin.
We also assume an elastic medium. Volcanic regions are known to host inelastic processes such as seismicity, damage, thermoplasticity, infiltration of and alteration by hydrothermal and magmatic fluids, that can affect both stresses and dike propagation. In particular, these inelastic processes compete with stress-generating mechanisms by homogenizing stresses (e.g., McGarr & Gay, 1978;Savage et al., 1992;Stephansson, 1988). Repeating magmatic intrusions may also bring the state of stress to isotropic in the long run: since they tend to open perpendicularly to 3 , the strain they cause tends to bring σ 3 closer to σ 1 (Bagnardi et al., 2013;Chadwick & Dieterich, 1995;Corbi et al., 2015Corbi et al., , 2016. Additionally, faulting and earthquakes may dissipate shear stresses over time. In other words, the stress contributions in Equation 1 can change or be altered. An accurate calibration of the stress state needs to take into account the relaxation of each stress contribution over time and space, discriminating between stress sources (in particular topography-altering events) that became active at different times. These processes are difficult to constrain and are currently accounted for through rough approximations. For instance, some works set the deviatoric stresses arising from gravitational loading of the edifice to zero Heimisson et al., 2015). Corbi et al. (2015) found that superposing the effect of caldera unloading to a volcanic edifice where the state of stress is set to isotropic, rather than fully loaded, better explained the orientation of eruptive fissures at Fernandina, Galápagos. Here we neglected such processes by creating scenarios where dikes propagate below and around a caldera but not within an edifice, as the height of all topographic highs in our scenarios (Section 2.3) is lower than or comparable to the MDT (Section 2.2.1).
As shown in Figure 1c, stress contributions of magma reservoirs are dominant only in the proximity of the stress source. Such effect, nonetheless, can be important in determining nucleation points for dikes (Grosfils et al., 2015;A. Gudmundsson, 2006), that we do not model precisely here, as well as attracting or repelling incoming dikes if the reservoir pressure is increasing or decreasing, respectively (Pansino & Taisne, 2019).
Stress contributions due to previous large earthquakes may also deviate dikes or arrest their propagation. This has been considered both through theoretical (Maccaferri et al., 2014 and analog (Le Corvec et al., 2013) modeling. The fault-generated stresses do not influence dike trajectories significantly unless they come to close proximity (e.g., Maccaferri et al., 2014). However, Maccaferri et al. (2016) showed how an incoming dike can trigger the slipping of a pre-stressed fault, and be stopped by the resulting compressive stress. Faults and dikes may also interact with each other, for instance alternately accommodating tectonic extension (Gómez-Vasconcelos et al., 2020).
All these stress sources can be integrated in our models as they stand now. Including stress mechanisms that are not well-constrained, however, ultimately adds more uncertainty to a model rather than improve it.
One major simplification in SAM is that of linear pressure gradients across the plane of SAM cracks (Section 2.2.1). SAM, as a simplified model, cannot deal with stresses that are too heterogeneous, although in the example shown in Lateral-Dike (Figures 3b-3d) it well compared to TIM, which can deal with stress heterogeneity at the scale of the individual triangular dislocations composing the dike meshes. An additional issue, not discussed here, is the potential heterogeneity in the dike internal pressure arising from the viscous flow of magma (Lister & Kerr, 1991) or pockets of bubble-rich magma within the dike (Costa et al., 2009). Non-linear gradients in the internal pressure may affect the direction of propagation of SAM dikes. In this regard, the analytical model by Pollard and Townsend (2018) computes the stress intensity factor at the tip of a 2D vertical crack under arbitrary distributions of normal tractions, and may be used in future works to estimate the error in K when using the linear pressure gradient assumption in SAM.
The outcomes of our synthetic scenarios show that dikes are deflected away from topographic lows (calderas), and attracted by topographic highs (hills, resurgent domes), even small-sized ones (e.g., the resurgent dome in Elliptic-Caldera). This is consistent with previous dike propagation and stress models considering topographic loading/unloading (Corbi et (Gaete et al., 2019;Mantiloni et al., 2021). The few synthetic scenarios we present here, however, are not designed to reproduce the wide variety of vent patterns observed at real calderas. They do, nonetheless, reproduce some common features of vent distribution in calderas. When a coastline is involved in our scenarios (Figures 5b-5e), most or all dikes end up on the mainland. This is compatible with vent patterns in similar natural settings, such as Campi Flegrei (Smith et al., 2011) or Aira caldera, Japan (Geshi et al., 2020). In our tests, no dike trajectories end up within the caldera, except in Refilling-Caldera and Elliptic-Caldera. Cases where past eruptive vents lie predominantly at or outside the caldera rim include most Galápagos volcanoes (Chadwick & Howard, 1991) and Aira caldera, Japan, (Geshi et al., 2020). Vents opening within a caldera can be observed in several other settings, like Newberry caldera, Oregon (MacLeod et al., 1982), Santorini caldera, Greece (Sigurdsson et al., 2006), or Campi Flegrei caldera, Italy (Smith et al., 2011). Intracaldera vent openings are predicted when the caldera is very shallow, unloading is reduced by refilling (Refilling-Caldera), or a resurgent dome is present (Elliptic-Caldera). Nonetheless, these factors are not always associated with intracaldera vents in nature (e.g., no eruptions have occurred at Long Valley caldera's resurgent dome after doming inception, Hildreth, 2004). Applying a model to a real caldera entails a deeper understanding of its evolution, stratigraphy and eruptive history, and requires dedicated work. For this reason, we chose not to apply our models to real calderas in this work, as running our model for a real scenario without a proper calibration of the stress state is no different than setting up a synthetic scenario with arbitrary stress. The fast dike propagation model we presented here is particularly suited for stress calibration procedures, such as the one by Rivalta et al. (2019). This will be the subject of future work.
Our model does not consider the viscous flow of magma within dikes and, as such, does not model dike velocity. The two approaches may be integrated by combining the pathways predicted by our model with existing models of dike velocity Pinel et al., 2017) or growth, such as Zia and Lecampion (2020), introducing a numerical model of propagation of planar 3D hydraulic fractures, or Möri and Lecampion (2022); Pansino et al. (2022). We also remark that different magma compositions may involve large differences in magma viscosity and density, and neglecting the viscous flow may undermine the predictive power of our dike models in case of high-viscosity magmas.
In both SAM and TIM, dikes are assumed to break away from the magma reservoir after nucleation, as dike propagation is entirely driven by external stress and magma buoyancy force. In a more general case, the dike may be coupled to a reservoir, as past dike intrusion episodes have suggested (e.g., M. T. Gudmundsson et al., 2016;Maccaferri et al., 2016). The direction of dike propagation, however, may still be controlled by the gradient of internal pressure and external stress rather than the pressure imparted by the chamber, even though accounting for the viscous flow may change that. Analytical models of propagating dikes coupled with a magma chamber (Segall et al., 2001;Rivalta, 2010;M. Townsend & Huber, 2020) are only available for fixed dike orientations and, as such, cannot be applied to 3D dike trajectories. In our context, increasing the volume of a TIM dike as it advances could be a rough approximation of a dike-magma chamber coupling. Our results from comparing TIM and SAM (Section 2.4) suggest that the trajectories would not differ much even for large volumes of TIM dikes, as long as they start aligned to the external stress field (see Figures 3a, 4a, and 4c). Including dike-reservoir coupling in SAM or TIM, however, requires dedicated work.
In conclusion, we have developed a fast and flexible dike propagation model, complementing the numerical model by Davis et al. (2020Davis et al. ( , 2021 Stress models, however, are still critical and not yet fully understood. In a real-case application, our scenarios would be the end point of a stress calibration, whereby the stress state of a volcanic region is constrained through a statistical procedure aiming at matching dike simulations with observations, such as past vent locations , orientation of exposed dikes (Maerten et al., 2022) or focal mechanisms (Zhan et al., 2022). Our model is well-suited for such purpose. Once the stress is calibrated, it may be used to perform a long-term forecast on future vent locations, while the more sophisticated model may be employed to produce short-term propagation scenarios for incipient dike intrusions.