The magmatic system under Hunga volcano before and after the 15 January 2022 eruption

One of the largest explosive eruptions instrumentally recorded occurred at Hunga volcano on 15 January 2022. The magma plumbing system under this volcano is unexplored because of inherent difficulties caused by its submarine setting. We use marine gravity data derived from satellite altimetry combined with multibeam bathymetry to model the architecture and dynamics of the magmatic system before and after the January 2022 eruption. We provide geophysical evidence for substantial high–melt content magma accumulation in three reservoirs at shallow depths (2 to 10 kilometers) under the volcano. We estimate that less than ~30% of the existing magma was evacuated by the main eruptive phases, enough to trigger caldera collapse. The eruption and caldera collapse reorganized magma storage, resulting in an increased connectivity between the two spatially distinct reservoirs. Modeling global satellite altimetry–derived gravity data at undersea volcanoes offer a promising reconnaissance tool to probe the subsurface for eruptible magma.

) with Dr = 0.5 g/cm 3 , for different combinations of inversion norms (smoothest, A, to most compact, I) using SimPEG.The inversions with the smoothest norms (A to C) attribute a very small density contrast to a large area and fail to recover the prism depth, shape, and density contrast, therefore we do not consider them for our data inversions.Preferred norms used in modeling Hunga volcano data are D, E, and G (Fig. S7 to S13).S4) with Dr = 0.5 g/cm 3 , for different combinations of inversion norms (smoothest, A, to most compact, I) using SimPEG.Source depth, shape, and density contrast recovery are more accurate for VGG inversions (Fig. S4) therefore we use VGG data to model Hunga volcano gravity anomalies.1,2,2,2], [1,1,1,1],  and [0,2,2,2

Fig. S1 .
Fig. S1.Timeline of marine gravity anomaly.Annual updates of the marine gravity grid (colorbar) show that no significant change is observed until after the 2022 eruption.Histograms of the yearly difference and uncertainty (pink) estimated for the difference in gravity anomaly.Eruptive episodes shown as red bar.Green arrow on last histogram shows the significant gravity decrease observed between August 2021 and August 2022.

Fig. S2 .
Fig. S2.Choice of crustal correction density for the Bouguer anomaly calculation.(A)Relationship between residual Bouguer anomaly and bathymetry for the optimal correction density of 2750 kg/m 3 and best fit least-square regression (orange line).(B) Least-square slope estimates (L2 slope), and (C) Pearson correlation coefficient (r) for correction density varying from 2200 to 2900 kg/m 3 .The correlation is minimized (closest to 0, shown as green line) for a correction density of 2750 kg/m 3 .

Fig. S3 .
Fig. S3.Synthetic block model for the VGG inversion.(A) Model mesh for a 4 km by 4 km prism buried at 6 km depth with Dr = 0.5 g/cm 3 .(B) Synthetic vertical gravity gradient (VGG) data due to the density anomaly in (A) with the addition of gaussian noise (1% of the maximum anomaly amplitude).

Fig. S4 .
Fig. S4.Density contrast recovery for a VGG synthetic model of a 4 km by 4 km prism at 6 km depth (Fig.S3) with Dr = 0.5 g/cm 3 , for different combinations of inversion norms (smoothest, A, to most compact, I) using SimPEG.The inversions with the smoothest norms (A to C) attribute a very small density contrast to a large area and fail to recover the prism depth, shape, and density contrast, therefore we do not consider them for our data inversions.Preferred norms used in modeling Hunga volcano data are D, E, and G (Fig.S7to S13).

Fig. S5 .
Fig. S5.Synthetic block model for the gravity (gz) anomaly inversion.(A) Model mesh for a 4 km by 4 km prism at 6 km depth with Dr = 0.5 g/cm 3 .(B) Synthetic gravity data due to the density anomaly in (A) with the addition of gaussian noise (1% of the maximum anomaly amplitude).

Fig. S6 .
Fig. S6.Density contrast recovery for the Bouguer gravity anomaly synthetic model of a 4 km by 4 km prism at 6 km depth (Fig.S4) with Dr = 0.5 g/cm 3 , for different combinations of inversion norms (smoothest, A, to most compact, I) using SimPEG.Source depth, shape, and density contrast recovery are more accurate for VGG inversions (Fig.S4) therefore we use VGG data to model Hunga volcano gravity anomalies.

Figure S11 .
Figure S11.Interpretation of gravity changes.Possible physical processes explaining the mass/density increase or decrease depending on the sign of the observed gravity anomaly before and after the eruption.Blue -/--indicates weak/strong negative anomaly, Red +/++ indicates weak/strong positive anomaly.0 indicates no gravity anomaly.

Fig. S12 .
Fig. S12.Depth slices through the 3D density model of the post-eruptive magmatic system.Inversion results with three different norm combinations: [1,2,2,2] (A), [1,1,1,1] (B), and [0,2,2,2] (C).Color scale shows the density contrast with respect to a reference density of 2.75 g/cm 3 .Caldera outline (dashed green line) and historical vents (stars).Thin black lines are the same as in Fig. S6 (shown for reference only), and dotted lines are model cross sections displayed on Fig. S13 and Figures 4 and 5 in the main manuscript.

Fig. S15 .
Fig. S15.Post-eruptive gravity data and uncertainty.(A) 1-arcminute marine gravity grid from satellite altimetry (V32.1)(21, 23) for the post-eruptive time interval (colored circles) superimposed on shaded relief of the pre-eruptive multibeam bathymetry.(B) Uncertainty on the marine gravity anomaly varies from 0.6 to 2.2 mGal with a mean of 1.3 mGal over the study area.
Fig. S18.Downsampling of bathymetry grids.(A) Pre-eruptive bathymetry grid with original 3-arcsecond resolution, (B) downsampled to 1-arcminute resolution using a Gaussian filter (standard deviation of 2 km), (C) post-eruptive bathymetry grid, and (D) downsampled to 1arcminute resolution used for the calculation of the gravity effect of the bathymetry interface to derive the Bouguer anomaly.