Coulomb failure of Earth’s brittle crust controls growth, emplacement and shapes of igneous sills, saucer-shaped sills and laccoliths

Abstract Tabular intrusions are common features in the Earth's brittle crust. They exhibit a broad variety of shapes, ranging from thin sheet intrusions (sills, saucer-shaped sills, cone sheets), to more massive intrusions (domed and punched laccoliths, stocks). Such a diversity of intrusion shapes reflects different emplacement mechanisms caused by contrasting host rock and magma rheologies. Most current models of tabular intrusion emplacement assume that the host rock behaves purely elastically, whereas numerous observations show that shear failure plays a major role. In this study, we investigate the effects of the host rock's Coulomb properties on magma emplacement by integrating (1) laboratory models using dry Coulomb granular model hosts of variable strength (cohesion) and (2) limit analysis numerical models. Our results show that both sheet and massive tabular intrusions initiate as a sill, which triggers shear failure of its overburden along an inclined shear damage zone at a critical sill radius, which depends on the emplacement depth and the overburden's cohesion. Two scenarios are then possible: (1) if the cohesion of the overburden is significant, opening of a planar fracture along the precursory weakened shear damage zones to accommodate magma flow, leads to the formation of inclined sheets, or (2) if the cohesion of the overburden is negligible, the sill inflates and lifts up the overburden, which is dissected by several faults that control the growth of a massive intrusion. Finally, we derive a theoretical scaling that predicts the thickness-to-radius aspect ratios of the laboratory sheet intrusions. This theoretical prediction shows how sheet intrusion morphologies are controlled by a mechanical equilibrium between the flowing viscous magma and Coulomb shear failure of the overburden. Our study suggests that the emplacement of sheet and massive tabular intrusions are parts of the same mechanical regime, in which the Coulomb behavior of the Earth's brittle crust plays an essential role.


Introduction
Subvolcanic processes in the shallow Earth's crust (e.g., magma storage, magma transport) govern the link between the plutonic and volcanic realms (Breitkreuz and Rocchi, 2018;Burchardt, 2018). Tabular igneous intrusions represent fundamental elements of igneous plumbing systems that accommodate the bulk of magma transport and emplacement in the Earth's brittle crust (Magee et al., 2017;Cruden et al., 2018;Galland et al., 2018). Many of these intrusions form in sedimentary basins, in which they may exert significant mechanical and thermal effects on hydrocarbon systems (e.g., Senger et al., 2017;Spacapan et al., 2018).
Tabular intrusions exhibit diverse shapes associated with contrasting thickness-to-diameter aspect ratios, ranging from thin sheets that have low aspect ratios (dykes, sills, cone sheets) (Kuenen, 1937;Hutton, 2009;Magee et al., 2017;Schmiedel et al., 2017b;Cruden et al., 2018;Galland et al., 2018), to thick, massive intrusions that have much higher aspect ratios (laccoliths, plutons, plugs) (Gilbert, 1877;Corry, 1988;Schmiedel et al., 2015;Mattsson et al., 2018). Presently, most theoretical and numerical models of tabular intrusion emplacement consider the intrusions' overburden as a bending thin elastic plate (e.g., Pollard and Johnson, 1973;Michaut, 2011;Galland and Scheibert, 2013). The mathematical formulation of these models builds on the following main assumptions: (1) the intrusion diameter should be at least 10 times larger than the   Galland et al., 2018). (b) Saucer-shapedsill with inclined sheet associated with shear failure at the intrusion tips (drawing from a field photo, Hutton, 2009). (c) Cone sheets, emplaced due to shear failure of the host rock (schematic modified after Kuenen, 1937). (d) Laccolith, magma emplacement is accommodated by shear failure where the strength of the host rock was exceeded. Shaded colors represent reconstructed eroded overburden (schematic drawing modified after de Saint-Blanquat et al., 2006). (e) "Punched" laccolith accommodated by shear failure of the overburden, magma body pierces through the overburden (profile modified after Breitkreuz and Mock, 2004;. (Red -magmatic intrusion, grey -weak sedimentary host rock, e.g. shale, yellow/orange -sedimentary host rock, e.g. sandstone, carbonate; For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) depth, (2) the thickness of the overburden should be 5 times larger than the thickness of the intrusion (see Ventsel and Krauthammer, 2001), and (3) the intrusion's overburden only deforms by elastic bending. Combining the first two geometric assumptions leads to the formulation that the intrusion diameter should be 50 times larger than the intrusion thickness.
These assumptions, however, are rarely fulfilled for laccoliths and felsic sills. Even if large mafic sills fulfill these assumptions, the geometric scaling of Cruden et al. (2018) shows that 90% of laccoliths and felsic sills are too thick (sometimes by one order of magnitude and more), to be addressed by elastic plate models. In addition, numerous field observations indicate that inelastic deformation (brittle faulting, ductile flow) accommodates for a significant part to the emplacement and growth of sills and laccoliths ( Fig. 1) (e.g., de Saint-Blanquat et al., 2006;Hutton, 2009;Wilson et al., 2016;Spacapan et al., 2017). Finally, Got et al. (2017) recently suggested that the effect of non-linear, inelastic damage is necessary to explain geophysical and geodetic measurements at Grimsvötn volcano, Iceland. The mismatch between geological and geophysical observations in the host rock around most sills and laccoliths and the assumption of elastic behavior questions the systematic application of these models to many tabular intrusions in nature.
Recent laboratory Montanari et al., 2017;Schmiedel et al., 2017a), numerical (Haug et al., 2017;Gerbault et al., 2018;Haug et al., 2018) and theoretical (Scheibert et al., 2017) models show instead that the inelastic brittle behavior of the Earth's crust might play a major role in the emplacement, growth and shape of both sheet and massive tabular intrusions. However, when these models are considered separately, they are not capable of predicting the extent to which elastic versus inelastic deformation are the dominant deformation mechanisms that accommodate the emplacement of tabular igneous intrusions. To overcome this challenge, we integrate quantitative laboratory experiments of magma emplacement in the brittle crust with limit analysis numerical modeling. Our laboratory experiments reproduce the geometry of several types of tabular igneous intrusions observed in nature, from thin saucer-shaped sills, cone sheets, domed laccoliths to punched laccoliths. Limit analysis numerical models reveal that the mechanism that initiates the emplacement of the laboratory intrusions is dominated by Coulomb (shear) failure of their overburden. While the final shapes and extents of thin saucer-shaped sills and cone sheets in our laboratory experiments can be explained by the development of a single inward dipping damage zone cutting the overburden, the shapes of thicker bodies, such as (punched) laccoliths, result from a more complex and pronounced deformation of the overburden. Our integrated approach shows how the inelastic brittle properties of the Earth's crust play a significant role in the emplacement of igneous intrusions in the brittle crust.

Laboratory experiments
The experimental method used in this study is similar to that presented by Schmiedel et al. (2017a). Our experiments aim at simulating basin-scale processes, thus one centimeter in our experiments represents between 100 m and 1 km in nature. The model rock and model magma consist of Mohr-Coulomb dry granular materials of variable cohesion C (Galland et al., 2006;Abdelmalak et al., 2016), and molten vegetable oil, respectively (Fig. 2). The granular materials used in this study consist of a mix of (1) a high cohesion fine-grained crystalline silica flour, which fails both in tension and shear (SF; C ∼ 600 Pa), and (2) low cohesion micro glass beads (GB; C ∼ 100 Pa; Fig. 2A). When mixed, they produce homogeneous mixtures (SF/GB) of intermediate cohesions that correlate linearly with the mixing proportions (Abdelmalak et al., 2016). These granular mixtures thus simulate a wide range of brittle rocks, from strong, competent rocks (e.g. limestone) to weak rocks (e.g. shale) ( Fig. 2A). The molten vegetable oil has a viscosity of η = 2 × 10 −2 Pa s at 50 • C (Galland et al., 2006). Guldstrand et al. (2017) and Schmiedel et al. (2017a) discussed that the oil is suitable for simulating magmas of intermediate to high viscosity in nature. The relevance of these materials and the scaling of the models are described in detail by Galland et al. (2014) and Schmiedel et al. (2017a) and references therein (see also Supplementary materials for further details). Fig. 2B shows a sketch of the experimental setup, which was first introduced by Galland et al. (2009) and modified in Schmiedel et al. (2017a). We prepared the experiments by compacting two layers of the same granular material in a 0.4 × 0.4 m box ( Fig. 2B; Table 1) using a high frequency compressed-air shaker (Houston Vibrator model GT-25). The compaction procedure was conducted as follows: (1) we poured a first layer of granular material of known mass in the box and gently compacted it, (2) we placed a flexible net on top of the compacted lower layer to simulate a thin horizontal heterogeneity, (3) we poured a second mass of the same Coulomb granular material and we compacted the whole model. We controlled the amount of compaction such that the net was located just in contact with the injection inlet and the density and so the thickness of the upper layer was known with great precision.  We estimate the density difference between the two layers of compacted granular materials to be <5%. This compaction procedure allows for a good reproducibility of the granular material density, and therefore a good control of the cohesion (Galland et al., 2009;Schmiedel et al., 2017a). We assume the flexible net in our laboratory models is a valid approximation to simulate a weak interface between sedimentary layers in nature, given that the strength contrasts in layered rocks can be as high as several orders of magnitude (Hoek et al., 1998;Schellart, 2000). However, the mechanical effect of the flexible net is hard to quantify. Note that the oil and the host granular materials exhibit distinct density contrasts (see supplementary material), which could po-tentially produce distinct buoyancy effects. However, Galland et al. (2009) showed that the overpressure contribution due to buoyancy was negligible with respect to overpressure driving oil propagation in the models. We hence infer that the buoyancy effects are negligible for interpreting our laboratory results.
In all our experiments we injected the oil with a constant volumetric flow rate Q of 40 ml min −1 (=5.1 × 10 −7 m 3 s −1 ). The diameter of the injection inlet d is 0.005 m, thus the average injection velocity of the oil is v i = 0.026 m s −1 . Because the net is located just at the inlet (Fig. 2B), the oil flows horizontally along the flexible net before propagating upwards. Thus, our experiments do not consider the feeder-to-sill transition. The horizontal extent of the flexible net was significantly larger than the final intrusion diameters, therefore boundary effects on the oil propagation due to the flexible net are negligible. The end of the experiments corresponded to the eruption of the oil; the durations of the experiments ranged between 4 and 88 s.
During the experiments, we monitored the surface deformation induced by the intruding oil using a photogrammetry system based on four synchronized cameras ( Fig. 2B) Schmiedel et al., 2017a). This photogrammetry system captured the model surface with a frequency of 1 Hz and produced highresolution (<0.1 mm) and high-precision (∼0.05 mm) topographic maps Schmiedel et al., 2017a). After each experiment, the oil solidified, we carefully excavated the intrusions in situ and we digitalized their 3D shapes using photogrammetry Schmiedel et al., 2017a). This procedure allows a direct comparison between the intrusion geometry and the associated surface deformation. Visual observations show some percolation of the oil into the granular host material during cooling. However, the percolation and the cooling take significantly longer than the actual duration of the experiments itself (Galland et al., 2006). Thus, we infer that the effects of oil percolation and cooling do not affect oil emplacement during the experiments.
Here, we present a series of laboratory experiments in which we vary two parameters: (1) the cohesion C of the model rock ranging between ∼100 and ∼600 Pa, and (2) the thickness of the overburden H ranging between 0.01 and 0.06 m (Table 1). Haug et al. (2017) propose that the transitions between initial sill (inner sill) and inclined sheets (outer sill) of saucer-shaped intrusions are controlled by shear failure of the overburden. The similarity between the saucer-shaped intrusions produced in the experiments of Galland et al. (2009), Schmiedel et al. (2017a, and in this study and the numerical models of Haug et al. (2017) leads to the following hypothesis: the transition between the horizontal initial sill and the inclined sheet occurs when shear failure of the brittle overburden occurs. However, in our experiments, we monitor the surface deformation only, as the granular materials are opaque and do not allow us to map subsurface structures associated with the emplacement of the oil. Thus, to test this hypothesis, we directly compare our experimental results with numerical models. We use the limit analysis software Optum G2 (Krabbenhøft and Lyamin, 2014). Limit analysis is a standard engineering tool commonly applied to assess the mechanical stability of slopes, rock masses and buildings (Sokolovski, 1960;Souloumiac et al., 2010;Haug et al., 2017). The principle of limit analysis is to predict approximate values of loads that bring a brittle solid to an imminent state of failure (Souloumiac et al., 2009(Souloumiac et al., , 2010Haug et al., 2017). Limit analysis is based on two theorems, the lower and upper bound theorems, also called collapse load theorems. The lower bound theorem is used to find the highest prescribed load where no failure occurs. The upper bound theorem is used to find the lowest prescribed load that ensures failure. This implies that the actual critical load leading to failure of a brittle solid is in between the loads calculated from the lower bound and upper bound theorems. If the range between load estimates calculated from the lower and upper bound is small, the actual critical load for failure is well constrained. The main output of the limit analysis simulations for this study is the total energy dissipation, an equivalent for the distribution of shear damage throughout the overburden (Fig. 2D).

Numerical simulations
In this study, we apply limit analysis to test whether the transitions between growing initial sills and the subsequent inclined sheets in our experiments can be explained by shear failure of the sill's brittle overburden. To test this hypothesis, we compare the laboratory inclined sheets with numerical shear damage zones induced by sills of the same diameters as the laboratory ones: if they match, we can conclude that weakening induced by shear damage and failure of the overburden controls the sill-to-inclined sheet transition and the growth of the inclined sheet. Thus, for each laboratory experiment presented in this study, we ran a limit analysis model using the same host properties and characteristics. The link between the experiments and the limit analysis simulations is straightforward: for each laboratory experiment, we take the initial mechanical conditions (density ρ, cohesion C and angle of internal friction φ of the granular material), the initial boundary condition (depth H ) and the radius of the base of the laboratory intrusion R, and we use them as initial and boundary conditions to run a limit analysis calculation. The series of limit analysis models thus match the experimental parameters of all laboratory experiments. Each limit analysis model considers a sill represented by an elongated pressurized cavity of radius R and depth H (Fig. 2C), whereas the sill thickness is set constant at 0.003 m. This value is an approximate value of the uplift measured above sill intrusions in the laboratory experiments of Schmiedel et al. (2017aSchmiedel et al. ( , 2017b, assuming that the uplift reflected the underlying sill thickness. Haug et al. (2017) show that the sill thickness is not a critical parameter as long as the thickness-to-radius (diameter) ratio is 1, which is generally the case for sill intrusions. The entire axisymmetric, numerical domain is 0.2 m wide and 0.1 m high similar to the laboratory box (40 × 40 × 10 cm). The system is solved using an adaptive grid with approximately 10 000 elements (Krabbenhøft and Lyamin, 2014). We use a free surface on the upper boundary, and a semi free-slip condition allowing only vertical displacement is applied to the right boundary ( Fig. 2C). At the lower boundary, a no-slip condition is used. Several tests have been run using different sill tip shapes and thicknesses, domain sizes and boundary conditions, which have all been documented to produce negligible effects (Haug et al., 2017.

Intrusion geometry
We observe variable intrusion geometries with respect to the experimental parameters. In the high cohesion experiments (100% SF, 50/50 GB/SF; Fig. 3; Table 1), the intrusions were typically saucer-shaped, with circular, flat, horizontal inner sill and outer inclined sheets (Fig. 3). In contrast, in low cohesion experiments ( Fig. 3; Table 1), the intrusions were massive, piston-like intrusions, with sub-vertical sidewalls (Fig. 3). Finally, in experiments with intermediate cohesions (80/20 GB/SF, 90/10 GB/SF; Fig. 3, Table 1), the intrusions exhibit mostly inclined sheets, with a subhorizontal bottom intrusion of radius smaller than those of intrusions emplaced in the high-cohesion experiments (<0.025 m; Table 1). The sheets in intermediate-cohesion experiments appear thicker than those in the high-cohesion experiments (Fig. 3). Note, however, that the thickness of the excavated intrusions includes the intrusion itself and the aureole of oil percolation, the thickness  (2) associated numerically modeled shear damage zone (black) calculated with limit analysis software OptumG2. Experiment 20 (Table 1) is a reproduction of experiment 12 with the same initial conditions and resulted in similar intrusion shapes, therefore only experiment 12 is shown here. of which cannot be quantified. In addition, oil eruption drained out the intrusion after the end of the experiment, therefore the final volume of the excavated intrusion is not the injected volume of oil before eruption. Hence, we will neither consider the intrusion thickness nor intrusion volume further. However, an estimate of the intrusion volumes prior to eruption was calculated using the experiment duration and the injection flow rate (Table 1). Fig. 3 shows significant variation of intrusion radius, ranging from 4.7 cm in Experiment 1 (highest cohesion; Table 1) to 0.8 cm in Experiment 4 (lowest cohesion; Table 1). In the log-log plot of Fig. 4A, the correlation between base intrusion radius-toemplacement depth ratio (R/H ) and cohesion-to-lithostatic stress ratio (C /ρ g H) exhibits a linear trend of slope 0.94, which can be approximated as a linear scaling. Note that the data from both sheet intrusions and massive intrusions plot on the same linear scaling law (Fig. 4A).
We also look at the correlation between the angle α which is either the average dip angle of the inclined sheets (sheet intrusions) or the dip angle of the sidewalls (massive intrusion), and the laboratory input parameters C /ρ g H and φ. Fig. 4B and 4C display that there is no clear correlation between the angle α and C /ρ g H or between α and the angle of internal friction φ observable. We note that the average values of α are higher for the massive intrusions (∼75 • ) than for the sheet intrusions (∼65 • ; Fig. 4B and 4C). In the models of Haug et al. (2018), the angle of internal friction φ and the cohesion C of the sill's overburden could be varied independently, so that it was possible to decipher the effects of each parameter on the angle α. In contrast, in our experiments, φ cannot be separated from C (Abdelmalak et al., 2016), so that the Figs. 4B and 4C are the results of the combined effects of φ and C .

Surface deformation
The surface deformation in the high-cohesion experiments (100% SF, 50/50 GB/SF; Table 1) shows a smooth dome structure with small h/w values (dome height h/dome diameter w) of ∼0.02 and small surface curvatures (Table 1). Additionally, orthorectified images of the model surfaces display small tensile fractures, which concentrate in the areas of the dome with the highest curvature (Fig. 5A) Fig. 5A; Table 1), the model surface deformation patterns exhibit domes, which combine characteristics of those of high-and low-cohesion experiments. The domes were not blocky (Fig. 6B), but the values of h/w and the curvatures of the domes were higher than those of high-cohesion experiments (Table 1).  The eruption of the oil, marking the end of each experiment, is always observed to be located along the outer edge of the surface deformation (e.g. Fig. 6B). Fig. 3 shows the comparison of the calculated damage zones and the characteristic intrusion profiles from our laboratory experiments. In experiments that produced sheet intrusions, the shapes of the calculated damage zones systematically match those of the inclined sheets (Fig. 3). In experiments that produced massive intrusions, the match between intrusion morphology and calculated damage shear zone is poor. For all experiments, the lowermost parts of the sidewalls appear sub-parallel to the calculated damage shear zone, whereas the upper parts are not. Note that in Experiment 3, the left sidewall matches well the calculated damage shear zone almost up to the top of the intrusion (Fig. 3). Note that the dip angles calculated from our simulations are systematically lower than those measured in our experiments, in particular for massive intrusions (Fig. 4).

Interpretation
Our experiments show that massive intrusions form in lowcohesion model rock. The high-amplitude surface deformation in these models exhibit local high curvatures (Figs. 3, 5B and 6C), which evidence dominant shear failure of the intrusions' overburden. Conversely, sheet intrusions form in high-to intermediatecohesion models (Figs. 3 and 6B). The low-amplitude and smooth surface deformation in these models evidence a different overburden deformation mechanism than in the low-cohesion experiments (Figs. 3, 5A and 6A). Tensile fractures observed at the surface suggest that at least part of the overburden's deformation is accommodated by tensile fracturing (Fig. 4A). Our results thus evidence the first-order effects of overburden brittle properties and behaviors on magma emplacement.
Both massive and sheet intrusions in our laboratory models exhibit a flat bottom (Figs. 2, 3 and 4). This suggests that they result from an initial spreading of the oil along the flexible net as a sill, as suggested for the emplacement of saucer-shaped sills (Malthe-Sørenssen et al., 2004;Galland et al., 2009;Galland and Scheibert, 2013) and laccoliths (e.g., Jackson and Pollard, 1990;Schmiedel et al., 2015). The key question in our experiments is: which mechanism(s) control the transition from this initial sill to either saucer-shaped sills or punched laccoliths? In low-cohesion experiments, the surface deformation associated with the emplacement of massive intrusions seems clearly related to shear failure, i.e. faulting (Figs. 5B and 6C), suggesting that the transition from the initial sill and the growth of the massive intrusion is controlled by shear failure of the overburden.
In contrast, surface deformation associated with the formation of inclined sheets in the high-to intermediate-cohesion models is not clear. The excellent match between the calculated damage zones and the laboratory inclined sheets (Figs. 3, 4B and 4C) strongly suggests that the emplacement of the inclined sheets is controlled by shear damage across the overburden. We thus infer that the transition between the initial sill and the inclined sheets is controlled by the shear failure of the overburden. These results suggest that the radius, and so the lateral extent, of both the sheet intrusions and the massive intrusions is controlled by the same process, i.e. the shear failure of the overburden triggered by the emplacement of the initial sill. This conclusion is supported by Fig. 4A, which shows that the scaled base radius R/H of both sheet and massive intrusions exhibit the same linear scaling with respect to C /ρ g H.
Another key question arising from our experiments is: what are the processes that lead to drastically distinct final intrusion morphologies? As shown above, the early mechanisms controlling the emplacement of the initial sill until shear failure of the overburden is common to both sheet and massive intrusions (Fig. 4A). Therefore the difference in terms of final intrusion shapes must result from later stage mechanisms controlling either the emplacement of inclined sheets or the inflation of a massive intrusion. In the sheet intrusion experiments, the low values of the aspect ratios h/w indicate that the uplifting intrusions' overburden is affected by limited strain. Finally, the orthorectified images of the deforming model surface exhibit only minor tensile fractures (Fig. 5A). We infer from these observations that the overburden dominantly fails along one main shear plane (Fig. 7), i.e. the reverse damage zone predicted by the limit analysis models (Fig. 3). The damage along this shear plane weakens the overburden, and controls the emplacement of the magma along inclined sheets, as suggested by Haug et al. (2017Haug et al. ( , 2018. In this scenario, the cohesion of the overburden is such that the uplifting block remains sub-rigid (except a few minor tensile cracks at the surface), which allows tensile opening along the early reverse shear planes to accommodate the emplacement of the inclined sheets (Fig. 7, right).
Conversely, in the massive intrusion experiments (low cohesion model rock), the walls of the intrusions do not match the slope of the damage zones calculated in the limit analysis models (Fig. 3).
One reason for such a mismatch is that limit analysis applies for small strains only and provides only a static solution for the initial distribution of shear damage when the overburden fails, but it does not simulate the evolution of the overburden deformation. However, the larger values of h/w show that the overburden is affected by large strains. In addition, the orthorectified images of the model surfaces show the formation of numerous, prominent fractures and faults that dissect the intrusions' overburden (Fig. 5B). The mismatch between the assumptions of limit analysis models and our observations suggests an additional mechanism of the overburden deformation. The structures visible on the orthorectified images resemble those observed in laboratory models of resurgent domes, which form as result of emplacement of massive intrusions (e.g., Acocella et al., 2000). In these models, two main sets of faults accommodate the uplift of the dome as a response to magma intrusions: (1) an early inward-dipping reverse fault, and (2) a subsequent outward-dipping normal fault, which accommodate the stretching of the overburden during uplift. In addition, these models show that the subvertical walls of the intrusions are controlled by the combined movements of the reverse and normal faults (e.g., Acocella et al., 2000). We infer from these resurgent dome models that both inward-dipping reverse faults and outward-dipping normal faults affect the overburden of the massive intrusions in our experiments, and that the subvertical walls of the massive intrusions are controlled by the combined activity of reverse and normal faults (Fig. 7, left). We argue that this complex faulting pattern occurs when the cohesion of the overburden is so low that subsurface tensile opening is not possible, in contrast to our sheet intrusion experiments.
To summarize, our laboratory and numerical results show that the initial emplacement of both sheet and massive intrusions is governed by the emplacement of an initial sill until reverse shear failure of the overburden occurs (Fig. 7). The final diameter of the sill is controlled by the onset of shear failure of the overburden. The subsequent emplacement follows two possible scenarios: (1) if the cohesion of the host is high enough, the magma is intruded along the weakened shear damage zones, which open in tension and lead to the formation of inclined sheets, or (2) if the cohesion of the host is low, the magma pushes upward its overburden, which develops secondary normal faults. The combined movements of the reverse and normal faults control the growth of piston-shaped intrusions, with steep walls that are not parallel to the early reverse damage zone (Fig. 7).

Discussion
Our results show that Coulomb shear failure of the overburden dominantly controls the emplacement and evolution of massive laccolithic plugs (Figs. 5B and 6C). This corroborates well with field observations (Breitkreuz and Mock, 2004;Schmiedel et al., 2015;Mattsson et al., 2018) and recent laboratory models (Montanari et al., 2017;Schmiedel et al., 2017a). The occurrence of reverse and normal faults to accommodate the emplacement of massive intrusions in the experiments is in good agreement with the resurgent dome models from Acocella et al. (2000). Note that the initial reverse faults in these experiments (Acocella et al., 2000) correspond to the peripheral reverse fault in our laboratory experiments and the shear damage zone in our numerical simulations (Figs. 3, 5 and 7). Our results thus confirm the conclusions of Schmiedel et al. (2017a), who argue that the established models of laccolith emplacement based on pure elastic bending of the overburden (e.g., Michaut, 2011;Currier et al., 2017) might apply only to a minority of laccolithic intrusions in nature.
In addition, our models show that shear failure also controls the formation of inclined sheets of saucer-shaped and V-shaped intrusions. These results corroborate the field observations that document significant inelastic deformation accommodating sill propagation (e.g., Kuenen, 1937;de Saint-Blanquat et al., 2006;Hutton, 2009;Galland et al., 2018) seismic data (Magee et al., 2017;Schmiedel et al., 2017b) and the numerical models of Haug et al. (2017Haug et al. ( , 2018. However, our models contradict established models for saucer-shaped sill emplacement (e.g., Malthe-Sørenssen et al., 2004;Galland and Scheibert, 2013). In these latter models, tensile stresses due to the elastic bending of the overburden deflect the tip of the propagating initial sill leading to inclined sheets. We propose instead that the propagating initial sill triggers shear failure of the overburden; the resulting inclined damage zone is a weakness used by the subsequently intruding magma (Fig. 7). Note that the conclusions of this study differ from those of Schmiedel et al. (2017a), who concluded that the emplacement of the saucershaped intrusions dominantly resulted from elastic bending of the overburden.
Our models show that the brittle shear failure of the overburden controls the diameters of both saucer-shaped intrusions and laccoliths through the common scaling law of Fig. 4A. Thus, we provide an alternative mechanism that differs from models proposed by Chanceaux and Menand (2016), Currier and Marsh (2015) and Currier et al. (2017), where the lateral extent of intrusions is controlled by the syn-emplacement cooling effects of the magma. Note, however, that the models of these authors account for purely elastic host rock (i.e. gelatin), whereas we argue that the assumption of a purely elastic host is likely not realistic for the emplacement of laccoliths and plutons. We therefore argue that shear failure of the overburden is a simple and realistic mechanism that controls the lateral extent of both igneous sheet and massive intrusions. Fig. 4A shows that the dimensionless diameter R/H of both saucer-shaped and laccolithic intrusions simply scales linearly with C /ρ g H. This scaling seems to differ from that proposed by Haug et al. (2018) (Equation (7)). However, a direct comparison is to be avoided, given that the boundary condition in the models of Haug et al. (2018) consider the magma pressure, whereas the boundary condition in our models consider constant injection flow rate, with variable magma pressure both spatially and temporally. In order to better integrate our laboratory models with the numerical models of Haug et al. (2018), measurements of the magma pressure during the experiments should be implemented.
Our laboratory experiments highlight the effect of the overburden's cohesion on intrusion morphology and emplacement, where low cohesion leads to massive intrusions and intermediate to high cohesion leads to sheet intrusions (Figs. 3, 4B and 4C). Nevertheless, field observations show that there is also a correlation between magma composition (i.e. viscosity) and intrusion morphology: (1) the majority of sheet intrusions are of mafic to intermediate composition, hence, low to intermediate viscosity magma, whereas (2) most of the massive intrusions are of intermediate to felsic composition with an intermediate to high viscosity , and references therein). The combination of our laboratory results and geological observations suggest that the development of igneous intrusions, and so their final shapes, is likely controlled by the balance between (1) the ability of a host rock to create open pathways for magma and (2) the ability of the magma to flow (Galland et al., 2014).
To test this hypothesis, we derive a simple analytical model to assess the mechanical behavior of the deforming brittle solid and the flowing viscous liquid, and how this controls the intrusion shapes in the experiments. The model considers the balance between (1) the mechanical work accommodated by the brittle deformation and uplift of the overburden W host and (2) the mechanical work accommodated by viscous dissipation in the flowing magma W magma . The supplementary material describes the details of the derivation of W host and W magma : where ρ is the density of the host and d the diameter of the inlet, respectively. The model assumes that the development and growth of intrusions result from a balance between the work of the viscous dissipation, and the work resulting from energy dissipation due to uplift and failure of the overburden. Such work balance implies deriving an expression of the intrusions' aspect ratio T /R as a function of the model parameters: This expression predicts that T /R scales with Π to the power of −1/3. In other words, it shows how low viscosity magma intruding into cohesive rocks is expected to produce thinner intrusions than high-viscosity magma intruding into weak rocks. This prediction is consistent with geological observations (Cruden et al., 2018 and references therein;Mattsson et al., 2018) and with existing models of sills and laccoliths emplacement (Roman-Berdiel et al., 1995;Currier and Marsh, 2015). We recall that during our experiments we cannot measure the thickness and radius of the intrusion, thus in the following analysis we use the aspect ratio of the uplifted dome h/w as a proxy of T /2R (see discussion by Schmiedel et al., 2017a).  (3). This match validates the physical relevance of our model to sheet intrusions, and quantifies how the morphology of tabular intrusions is controlled at first-order by a mechanical equilibrium between the viscous flow of the magma and the Coulomb brittle deformation of the overburden. In the population for the massive intrusions, the values of h/w plot significantly higher than the values defined for the sheet intrusions (Fig. 8). The fact that the massive intrusions do not agree with our theoretical prediction suggests that another more complex mechanism of overburden's deformation accommodates the emplacement of massive intrusions. This is in good agreement with our interpretation that several structures, i.e. reverse and normal faults, dissect the overburden and control the growth of massive intrusions. Another explanation is that the massive intrusions are so thick that the magma flow regime, which is assumed to fulfill the lubrication theory in our theoretical model (see supplementary material; Bunger and Cruden, 2011) is not valid because intrusions are significantly thicker.

Conclusions
This paper describes the integrated results of laboratory and numerical modeling of tabular intrusions emplaced in the Earth's brittle crust. Our laboratory models are able to simulate the emplacement of both sheet intrusions, i.e., saucer-shaped sill, cone sheet, and massive intrusions, i.e. punched laccoliths (Schmiedel et al., 2017a). In this study, we tested the effect of host rock's Coulomb brittle properties, and specifically the failure mode, on the emplacement of magma. The main conclusions from our study are the following: 1. Saucer-shaped sills, cone sheets and laccoliths all initiate from an initial sill. When the initial sill reaches a critical diameter, it triggers shear failure of its overburden along inclined inwarddipping reverse shear damage zones. 2. The subsequent evolution of the intrusion can follow two scenarios according to the cohesion of the overburden. In highcohesion models, the intrusions develop inclined sheets that follow the inward-dipping shear damage zones, which accommodate magma flow by dominant tensile opening. In lowcohesion models, the intrusions develop as massive pistonshaped bodies by pushing upward their overburden, which is dissected by both reverse and normal faults, inwards and outwards dipping. 3. The dimensionless diameter R/H of the initial sill, and so the diameter of the intrusion, is linearly correlated with the dimensionless cohesion C /ρ g H of the overburden (Fig. 4A).
4. We propose a simple analytical model that predicts the thickness-to-radius T /R aspect ratio of the intrusions as a function of our model parameters. This theoretical model considers the mechanical equilibrium between the viscous flow of the magma and the Coulomb brittle deformation of the overburden. The well predicted scaling of the thickness-to-radius T /R aspect ratios of the sheet intrusions in our laboratory experiments validate the physical relevance of our theoretical model for the emplacement of saucer-shaped intrusions and cone sheets. Conversely, our theoretical model does not yet apply to massive intrusions, the emplacement of which is accommodated by a more complex deformation mechanism of the overburden. 5. Our results suggest that the emplacement of sheet and massive tabular intrusions are parts of the same mechanical regime, in which the Coulomb behavior of the Earth's brittle crust plays an essential role.