Southward growth of Mauna Loa’s dike-like magma body driven by topographic stress

Space-geodetic observations of a new period of inflation at Mauna Loa volcano, Hawaii, recorded an influx of 0.11 km3 of new magma into it’s dike-like magma body during 2014–2020. The intrusion started after at least 4 years of decollement slip under the eastern flank creating > 0.15 MPa opening stresses in the rift zone favorable for magma intrusion. Volcanoes commonly respond to magma pressure increase with the injection of a dike, but Mauna Loa responded with lateral growth of its magma body in the direction of decreasing topographic stress. In 2017, deformation migrated back, and inflation continued at the pre-2015 location. Geodetic inversions reveal a 8 × 8.5, 10 × 3 and 9 × 4 km2 dike-like magma body during the 2014–2015, 2015–2018 and 2018–2020 periods, respectively, and an average decollement slip of ~ 23 cm/year along a 10 × 5 km2 fault. The evolution of the dike-like magma body including the reduction in vertical extent is consistent with a slowly ascending dike propagating laterally when encountering a stress barrier and freezing its tip when magma influx waned. Overall, the magma body widened about 4.5 m during 2002–2020.

northward migration of the inflation source rather than an instantaneous shift, consistent with the gradual solidification of magma in the tip of the dike.
Figs. S1.5 addresses the question whether the May 4 2018 Kilauea earthquake impacted the displacement pattern. The May 2018-May 2020 velocity (without the earthquake, Figs. S1.5i,j) is identical to the April 2018-May 2020 ( Fig. S1.4g,h), indicating that the influence of the earthquake is negligible.

Modelling approach
We model the magma sources as rectangular, uniform opening dislocations 48 and point sources 49 , and the decollement fault underneath the volcanic pile as a shear dislocation, all in an isotropic elastic half space (Poisson's ratio of 0.25). We use a Bayesian inversion method to determine the best fitting source models which minimize the residual between the observed data and model predictions. We sample the LOS velocity maps and reduce the number of data points using an adaptive quadtree approach 50 . For each dataset we compute a variance-covariance matrix using an experimental semi-variogram 51 which is used to weight the data points using the inverse of their variance. We only use the horizontal GPS velocities owing to high noise levels in the vertical. We use the Geodetic Bayesian Inversion Software (GBIS) 52 , which uses a Markov Chain Monte Carlo approach with Metropolis-Hastings algorithm. We account for the elevation effect of the topography using reference elevations obtained by comparing the displacement fields from a halfspace model with those from computationally expensive varying-source depth models 18 ; the reference elevations for each time period are given in Table S2.3).
To fit the data, we consider three sources: (i) inflation of a dike-like magma body in the rift zone, (ii) inflation of a magma chamber under the caldera, and (iii) slip along the basal decollement, represented by a uniform opening dislocation (referred to as dike), a Mogi source, and a uniform, slightly upward dipping dislocation, respectively, all embedded in a homogeneous elastic halfspace. We also tested the Compound Dislocation Model (CDM) 53 , but obtained poorer fits. We consider various model configurations consisting of one or multiple sources and invert for the location, geometry and intensity (opening displacement, volume change and dip slip) with some constraints. For the dike we assume that it is vertical (dip fixed at 89.5°). For the decollement slip we assume uniform slip in direction perpendicular to the dike. As the data don't have the resolution to constrain the fault's location and size, we assume a 10*5 km 2 -sized patch, constrained to be located along the basal decollement fault (at a depth of 10 km b.s.l. under the caldera, 5° upward dipping towards east, parallel to the dike). The best-fitting models are characterized by a minimum in Root Mean Square Error (RMSE) between data and model predictions, calculated as: where d is data, m is model displacement and w is the weight of GPS data to InSAR data (1.0 in our case) and N is the total number of observations. Assuming suitable uniform priors (Table S2.2), we used a Markov Chain Monte Carlo based Bayesian inversion algorithm GBIS to derive posterior probability density functions (PDFs) for the model parameters. From the PDFs, we derive the optimal, mean and variance (2 ) for each model parameter. We estimated the model parameters using 1 million iterations and discarded the first 50000 iterations as initial burn-in 52,54 .

Modelling results
We consider five different time periods, the three time periods discussed in the paper for which we have InSAR and GPS data ( For each time period, we sub-sampled the InSAR data to around ~400 samples for each dataset. For the semi-variogram 51 used to construct the variance covariance matrix, we used a sill and nugget of 1x10 -6 m 2 and 1x10 -8 m 2 respectively. All the sources are allowed to vary within a set of non-informative priors (see Table S2.2).
The model fits for time periods 1,2,3 and are shown in Figs. S2.1 -S2.3. The model fit for time period B is shown in Fig. 2v in the main paper. The best-fitting model parameters along with variance (95 % certainty) are summarized in Table S2.1. The marginal and joint probability density distributions for models B and C for time periods 1,2 and 3 are shown in Figs S2.6 to S2.11. For comparison between different configurations, we also report the potencies for dislocation (surface area*opening), decollement fault (surface area * slip) 55 and mogi source (1.8*volume change) 16 .
For the three periods, we can discard Model A because it is associated with LOS increase (subsidence) over the dike which is not observed (Fig. S2.1 -S2.3). The addition of a mogi source (Model B) leads to a better fit to the data. Our preferred configuration is model C which includes decollement slip. The dike is located under the summit caldera parallel to the rift and the mogi source is located under the eastern caldera rim. The potency of mogi source (see caption to table S2.1) is 40-80 % of the dike. Although models B and C have similar RMSE, we select model C because the addition of the decollement slip explains the GPS-observed east flank motion and is consistent with the elevated seismicity along the decollement (Fig. 1e).
The marginal and joint probability density distributions for the model configurations (models B and C) show that the model parameters are well constrained (Fig. S2.6-2.11) and in both configurations there is a negative correlation between dike width to dike depth and dike opening, as well as a positive correlation between dike opening and dike depth. Furthermore, for all three time periods model configuration C shows a positive correlation between the volume change of the mogi source and fault slip. 1D histograms of the location of the Mogi source and the dike vertical extent elucidate the changes between time periods (Fig. S2.12 and S2.13).
To compare our results with the previous intrusion (2002-2005, time period A), we modelled the data from Amelung et al. 10 . The model fit is presented in Fig. S2.4.
The geodetic data don't have the resolution to resolve the precise location of slip and/or temporal variations as the net velocities at the surface are less than 1 cm (see Fig. S2.5). The same surface displacements are produced for example by 1.5 meters slip along a 10*5 km 2 patch (as modelled) or by 0.25 m slip along a 30*10 km 2 patch (for more see table S3.5.1).

Reference elevation for the model depths
As we use elastic half space models, we have to account for the elevation of the data points which varies considerably. Williams and Wadge (1998) 18 do this by (i) adjusting for each data point the source depth based on a point's elevation and (ii) using the mean elevation of all data points as the reference elevation. This method is referred to as the depth-adjusted approach. Inversion using this method is computationally expensive as the Green's functions need to be calculated at every point at every iteration. We found that the obtained model parameters using this approach and the depth-not-adjusted approach are nearly identical, except for model parameters describing the source depths. We therefore used the depth-not-adjusted approach for model configuration A to estimate a reference elevation such that depth a.s.l. from both, the depth-not-adjusted and the depth-adjusted inversions are the same. This reference elevation is then used for all model configurations. We found for the three time periods reference elevations of 3325, 3939, and 3894 m above sea level respectively (Table S2.3).

Summary of source strengths
The strengths of the sources are described by the opening rate, slip rate, the lengths and widths of the dislocations, and the volume change rates of the Mogi source (Table S2.1). The corresponding potency rates for the five time periods are given in table S2.5. In order to obtain excess pressure rates on the complex magma body, we convert the potency rates into equivalent opening rates over the modelled dislocation, which then can be used to calculate the excess pressure rates (Table S2.4) (see section S3.1 for definitions). As this excess pressure is dependent on the size of the dislocation which varies over the three time periods, we also calculated excess pressure over an 9x5 km 2 dislocation. The two different configurations are referred to in the table as varying-crack size and uniform-crack size model, respectively.

S3.1: Finite element method
To calculate the stress field under Mauna Loa, we built a Finite Element Model (FEM) as shown in Fig. S3.1. The topography of the Big Island is built into the model using SRTM dem (90 m). Boundary conditions include roller supports on all sides and fixed boundary at the bottom (Fig. S3.1a) following 57,58 . The top surface is let free to deform in all directions. The mesh for the computation is fine near the sources of deformation and coarsens towards the outer boundaries (Fig. S3.1b). Total size of the model is 300 km x 300 km x 300 km consisting of >100000 triangular elements. The size of the model is big enough in all directions so that calculated displacements are not affected by the boundaries. The section for the calculations is taken on a line passing through the modelled dike on the rift zone as indicated in (Fig. S3.1c). As described in the main paper, we compute perturbations to the lithostatic stress state under Mauna Loa due to, (i) magma body pressurization, (ii) decollement slip and (iii) topographic load.
We modelled our magmatic sources as lens shaped cavities in an isotropic medium and its perturbation to the lithostatic stress (magma body widening) is applied to the cavity walls as constant excess pressure ( ), where, b is the widening of the dike, E is Young's modulus and L is half of the dip-dimension of the dike. ν is Poisson's ratio.
Since our best fitting model is a complex magma body described by a combination of an opening dislocation and point source, we used the sum of dike and mogi potencies to compute equivalent widening (b) on a 9 * 5 km 2 dislocation. The thickness of the lens shaped cavity is kept at 50 m, previous studies have shown that this choice won't have a significant influence on the final result 58 .
We modelled the decollement as a 10*5 km 2 shaped rectangular cavity and a slip is applied to the walls. Topographic load is applied as body load on top of a lithostatic crust corresponding to ⍴gh (⍴ is the density; g is the acceleration due to gravity; h is the elevation above 1.7 km a.s.l.).
Stress perturbations are computed on a section passing through the magma body along the rift zone as shown in Fig. 3.1c.

S3.3 Analytical and FEM models comparison:
To validate our FEM results, a comparison is made using analytical dislocations 48 for 1m opening over a dike (Fig. S3.3a) and 1m dip slip over a thrust fault (Fig. S3.3b).

S3.4: Stress concentrations arising around the magma body
In Fig. 4a, c we show the rift-perpendicular normal stress perturbations along a section passing through the rift zone due to decollement slip and the topographic load, calculated assuming an isotropic elastic material. However, if the magma body is considered a fluid-filled cavity that deforms in response to imposed stress, it leads to stress concentrations along the boundaries 60 . The normal stress perturbations with cavity present are shown in Fig. S3.4a,c. For the cavity we used lens-shaped dike with dimensions of dislocation from time period 1 and maximum thickness of 50 m at center, gradually reducing to 0 m at edges. The choice of the cavity thickness does not significantly affect the stresses produced 58 ). However, stress concentrations near at edges of the cavity are ~50% higher when compared to without cavity.

S3.5: Earthquake potential along decollement fault
To evaluate the effect of inflation of the dike-like magma body on decollement slip, we considered three different-sized faults (10*5 km 2 , 20*20 km 2 , 30*30 km 2 see Fig. S3.5), divided them into 1*1 km 2 patches, and used the finite element model to calculate the slip in response to the imparted magma body widening (Table S3.5.1). We calculate a total moment ( 0 ) in Nm and Moment magnitude 61 ( ) as: A summary of moment imparted over the decollement due to magma body widening since 2002, is given in table S3.5.2.

S3.6 Tensional stress due to magma body widening in the rift zone:
Tensional stress due to magma body widening from three time periods (Fig. S3.6) used in the modelling and cumulative (Fig.  4d) accumulated so far in the rift zone has been calculated. Tensional stress in the rift zone accumulates as the magma body widens and gets relieved in the zone where the magma body intrudes 29 .       for a 1m slip on 10 km x 5 km Mw is 5.9, (**) calculated using = 2 3 log( 0 ) − 6.03 for 0.3 m slip on 10*5 km 2 is 5.55. log( 0 ) − 6.03.