Highlights

  • Quantum gravity technology on satellites fuels innovative Earth observation strategies, driving technological advancements

  • Submarine volcanoes lurk in oceans, undetected hazards. Enhanced gravity missions can spot giants like Fani Maoré (7.7 Gt)

  • MOCAST+ simulations with quantum gradiometer and clocks can improve seamount mass detection over GRACE

1 Introduction

The eruptions of submarine volcanoes alter the bathymetry, but these changes remain hidden below the ocean surface if the eruptions are non-explosive, and the observation through gravity is the only alternative of detecting these changes from remote observations. This work aims to define the noise level requirements of an improved mission to detect marine volcanic eruptions. We create a review and associated database of marine volcanic eruptions, starting from the Hunga Tonga Hunga Ha’apai (HTHH) explosive eruptions of 2014/2015 and 2021/2022, the silent Fani Maoré eruption of 2018–2021, describing further worldwide submarine volcano eruptions as the El Hierro 2011–2012 and Axial seamount 2015, and discuss their possible gravity signal. Submarine volcanoes are particularly difficult to reach with shipborne observations, due to their remote and sparse locations. Complete and systematic coverage of on-site observations is impossible, making satellite gravity observations particularly relevant and prone to give a significant scientific contribution. We test whether a next generation Quantum Space Gravimetry (QSG) mission carrying a payload of gradiometers and clocks in multiple satellites could have the sensitivity to detect these volcano eruptions.

The noise characteristics of an exemplary mission carrying quantum technology have been assessed for the MOCAST+ (MOnitoring mass variations by Cold Atom Sensors and Time measures) mission study, supported by the Italian Space Agency (ASI) (Rossi et al. 2023; Migliaccio et al. 2023). This proposal is a successor of the MOCASS (Mass Observation with Cold Atom Sensors in Space) mission (Migliaccio et al. 2019; Reguzzoni et al. 2021; Pivetta et al. 2021), consisting in a constellation of multiple satellites with a payload of quantum technology, time measurement and cold atom gravity gradiometers. We assess the expected improvements of this QSG mission with respect to the present technology GRACE-FO satellite mission (Kvas et al. 2019). The simulated gravity fields and associated errors of MOCAST+ are available in terms of spherical harmonic expansions, for which reason we take the spherical harmonic expansion products of the GRACE-FO gravity field observations for comparison (Level 2 data product). It can be argued that for GRACE-FO, a local analysis above the seamount with the lower level products (e.g., Level 1) as the intersatellite distance or intersatellite acceleration between the two satellites could lead to a lowered local noise level. We agree but cannot use such analysis, because it would not allow to compare the noise level with the MOCAST+ , which is not measuring intersatellite distance, and the only noise estimates which are available are in spherical harmonics. Therefore, the sensitivity assessment of MOCAST+ and comparison to GRACE is related to the Level 2 data products. Furthermore, the intersatellite distance processing requires a skilled geodetic dedicated personnel.

In the previous geophysical sensitivity assessment of the MOCASS satellite mission (Migliaccio et al. 2019; Pivetta et al. 2021), the first version of a seamount mass change estimate was devised—it involved gravity forward modeling and analysis in the spatial domain. While fit for purpose, that strategy runs into some limiting aspects, chiefly the comparison with the gravity field errors expressed in terms of error degree variance spectra. We here use a localized spectral analysis of the geophysical signals, which we describe below in general terms. In addition to that, we carry out the forward modeling of masses in the spherical harmonics (SH) domain, which saves additional synthesis–analysis steps and allows easily synthesizing scaled versions of the eruption volumes. Pivetta et al. (2022) evaluated the sensitivity of MOCAST+ to lakes, glaciers and subsurface hydrology. The present study complements that analysis with realistic submarine volcanic eruptions, implementing a novel detectability analysis method tailored to seamounts’ peculiar signal patterns.

The anticipated gravity signal from a submarine volcano exhibits temporal variations, either abruptly, as in the case of an eruption, or gradually, owing to sustained deformation and magma ascent. The temporal characteristics of this geophysical quantity are crucial for observations since the error in the gravity signal obtained from a satellite mission diminishes with an extended observation period. As the number of observations from the satellite increases over a more extended period, the root-mean-square error of the gravity field derived from observables decreases. The volcano eruption has a prolonged impact, notwithstanding sudden, step-function-like mass changes. Consequently, these phenomena can be detected by integrating observations over a much greater time interval than the event itself, as the integration time can be flexibly chosen both before and after the event. The chosen integration time affects the time resolution at which the phenomenon can be observed, with a longer observation time implying a coarser time resolution. For a signal to be discernible in the spectral domain, it must surpass the spectral noise of the mission's gravity product, requiring the signal spectrum to be above the error spectrum, at least within a range of spherical harmonic (SH) degrees.

The signals generated by mass movement in the atmosphere, ocean and hydrology (AOH) create aliasing phenomena, which have been estimated to be the limiting factor of near real-time products of gravity satellite missions, as they create an additional noise that is larger than the instrumental noise (Abrykosov et al. 2022). AOH signals are strong in the seasonal and yearly time variation (Dobslaw et al. 2015, 2016). This noise is reduced integrating the observation over time, which leads to a nulling average in the case of periodic signals (Abrykosov et al. 2021). Sufficient attenuation allows the error in the gravity models to approach the instrument-only noise contributions. Presently, the full estimate of the recovery of a seamount signal in the presence of aliasing has not been done yet, our study being preparatory to such a study, that of a synthetic gravity signal from a simulation including both AOH signal and the target seamount signal. The user requirements to a planned mission are formulated such to gain a significant improvement in the observation of a given geophysical phenomenon (e.g., Pail et al. 2015). This can be measured in terms of “completeness” of the detected phenomena, represented by the percentage of retrieved events for a given size. An improvement in completeness means recovering the signal down to small geographical extents or masses and/or small mass changes in time (i.e., slow longtime trends). Defining a database of geophysical signals and the spectral characteristics of the static and time-variable gravity signal they generate, contributes in setting the observational requirements and guiding mission scenarios and technological choices on the instrumentation for the next planned gravity missions. The Mission Requirement Document (MRD), e.g., Haagmanns et al. (2020) is set out in the planning stage of the mission, and includes the different applications, possible services to be developed with the satellite data, and the time and spatial resolution at target and threshold level of the mission. With the decision of the NGGM/MAGIC ESA-NASA mission to be in orbit in the late 2020s to early 2030s, the input from the user-communities is welcome and timely, and will be used to populate the Mission User Requirements Document. Initiatives include a public User-Requirements questionnaire (https://www.asg.ed.tum.de/en/iapg/qsg4emt/) and discussion forums, such as a Splinter meeting held at the EGU General Assembly 2023 (SP48 https://meetingorganizer.copernicus.org/EGU23/sessionprogramme/5179).

In this study, after giving a review on a selection of seamount eruptions, we define the expected mass change at the seamount, the associated gravity footprint and finally the detectability assessment of the MOCAST+ mission. A brief discussion of the sensitivity to solid Earth signals had been covered in Migliaccio et al. (2023) and includes HTHH. The mass estimate for the HTHH referred to the 2014/2015 eruption, in which a new island was formed, with an exposed surface of 1.8 km2, and an estimated volume change of 0.5 km3, based on documentation of Colombier et al. (2018). One month after that paper was accepted, the great HTHH eruption occurred (January 15, 2022), which caused an over tenfold volume change. We here discuss the much greater change of 2022.

In general, the mass increase in the bathymetry, retrievable from repeated shipborne depth-sounding acquisitions is an upper limit of the mass change; the net mass change could be lower if the magma chamber is not or only partially refilled from below.

2 Sensitivity Analysis Methods

We here describe the methods required to model the gravity signals of the selected seamounts and the formalism to determine the fields in spectral domain. We document the steps with which we obtain the mass variations associated with the volcanic eruptions, needed for the calculation of the gravity field.

2.1 Gravity Change Modeling

Our strategy to estimate the gravity signal due to the mass change of the volcano generated by a submarine eruption collocates the erupted mass on the summit of the volcano. We define the erupted mass from publications documenting the unrest episodes, choose a characteristic surface area for the seamount and position the erupted volume on the top of the seamount. The thickness of the erupted volume is obtained from the volume and the area, if no details are found in the literature. The position of the synthetic mass is chosen by centering the mass on the summit of the seamount’s bathymetry. Since the bathymetry data are affected by outliers (spikes) and some grid data are not defined (NaN) a nearest-neighbor interpolation function and a simple median filter are applied, obtaining a smooth morphology that makes it easier to distribute the synthetic mass over the seamount. We calculate the gravity field and its spectrum using the analytical expression of the spherical harmonic expansion of a spherical cone (Sutton et al. 1991) following the formalism reported in Supplementary material S1.

2.2 Localized Spectral Analysis in the SH Domain

The signals of interest have regional spatial distributions, meaning that they are spatially limited. Most of their power is concentrated in a limited geographical extent, their distribution over the Earth surface being inhomogeneous. The gravity field footprint in contrast is much larger than the diameter of the seamount which is favorable for detecting the mass change from a satellite, as long as the observation accuracy is sufficient and under the assumption that ocean and atmosphere noise has been correctly modeled. In formal terms, the global functions of the gravity signals we are studying are always spatially non-stationary. On the other hand, the error estimates on the SH coefficients of the global geopotential retrieved by a satellite mission are stationary, for the most part (Han and Ditmar 2008). Under these conditions, regional analysis is necessary to provide a consistent detectability assessment of the geophysical signals. As Han and Ditmar (2008) pointed out, performing a direct comparison between satellite spectral error estimates and global power spectra of signals within a “confined geographical regime” provides an incorrect picture. As we show in the example of a synthetic model resembling the Mayotte volcano in Section 2.3, the global SH spectrum of a local signal has much lower power per degree compared to a globally distributed signal with comparable amplitudes in the area of existence of the signal.

To avoid under-estimating the signal-to-noise ratio between the degree variance spectra of our signals of interest and the error degree variance of a simulated gravity field solution, we employ a regional analysis technique: the spatio-spectral localization as defined by Wieczorek and Simons (2005, 2007). The spatio-spectral localization, operating in the spherical harmonics domain, is based on the principle of applying a windowing function h to a global signal \(f\left( \Omega \right)\) the effect of which is deconvolved in the spectrum through the localization processing. In Supplementary material S2, the method is described, and moreover, the reader is referred to Wieczorek and Simons (2005) for a detailed derivation.

The cap radius \(\theta_{0}\) of the windowing function is one relevant parameter in the spectral localization and is commonly chosen so as to cover the signal of interest satisfactorily (Han and Simons 2008). This reflects its well-established use in the analysis of data on the sphere, even in the more general case of non-zonal arbitrary windows, which can be tailored to the study area (e.g., Harig and Simons 2016) when analyzing signals with more complex spatial patterns.

As discussed in Supplementary material S2, the cap radius \(\theta_{0}\) controls the spectral window bandwidth \(L_{{\text{h}}}\), limiting the lowest degrees of the calculated spectrum. There is a trade space between choosing a smaller cap radius and the omission of the low- and high-end of the spectrum that it entails. A small window may seem beneficial in retrieving the signal generated by small body; however, this does not hold true when most of the energy of the signal lies in the lowest degrees. As we show in Fig. 1, this is the case even for a small seamount-like source. As a metric to describe this omission, we define the ratio of gravity signal of the seamount under analysis integrated over a cap of given radius \(\theta_{0}\), with respect to the global signal. This is expressed by the following ratio, which we compute for the change in disturbing potential \(T\left( {\theta ,\phi } \right)\) due to the seamount mass change:

$$\Lambda \left( {\theta_{0} } \right) = \frac{{\mathop \smallint \nolimits_{0}^{2\pi } \mathop \smallint \nolimits_{0}^{{\theta_{0} }} T\sin \theta d\theta d\phi }}{{\mathop \smallint \nolimits_{0}^{2\pi } \mathop \smallint \nolimits_{0}^{\pi } T\sin \theta d\theta d\phi }}$$
(1)
Fig. 1
figure 1

Normalized signal amplitude and \(\Lambda \left( {\theta_{0} } \right)\) ratio of cumulative (cap integrated) signal retrieved for a varying radial distance from the seamount source. A Amplitude of gravity potential, normalized at the maximum value, for different gentle-cut truncations of the SH coefficients (\(l_{\max }\) the half-point SH degree in the tapered truncation, for a 24-degree wide taper). B Ratio of the area-integrated signal retrieved inside a cap of increasing radius \(\theta_{0}\), respect to the globally integrated signal (Eq. 1). The values on the legend and the arrows denote the cap radius at which the threshold is reached, above which more than 97.5% of the global signal is encompassed by the cap

Full size image

Therefore, a criterion can be defined on the basis of a chosen ratio \({\Lambda }\), using a cap with \(\theta_{0}\) that satisfies it.

In Fig. 1, we show the ratio \({\Lambda }\left( {{\uptheta }_{0} } \right)\) for a range of cap radiuses. It must be noted that increasing the maximum spatial frequency at which a non-bandlimited signal is synthesized—so the maximum SH degree at which a SH expansion is truncated—increases both its maximum amplitude (at zero distance from the signal source) and the amount of signal encompassed at the same cap radius (i.e., the signal peak becomes higher and narrower). We thus show the effect of truncating the SH expansion at different maximum degrees, which we implement in a tapered fashion to avoid excessive ringing with a tapered “gentle cut” filter as in Barthelmes (2008). Each maximum degree shown in Fig. 1 corresponds to the 0.5 value of the tapered filter, with the taper having a bandwidth of ± 12 SH degrees (e.g., \(l_{\max } = 35\) implies that the taper starts at \(l_{1} = 23\) and reaches 0 at \(l_{2} = 47\)).

We calculate the gravity field for a prototype seamount as the Fani Maoré discussed below, of diameter 0.028° (ca. 3.2 km) and mass 7.7 Gt, requiring that 97.5% of the footprint be covered by the cap of radius \(\theta_{0}\). The seamount is 820 m high, with its base at 3400 m depth (below sea level). The gravity field is calculated at sea surface. For a range of SH degrees consistent with the resolution of the gravity models provided by the mission configuration under test, values of \(\theta_{0}\) allowing for a \(\Lambda\) of 97.5% are in a range between 6.0 and 6.9 degrees. This shows that the footprint of the gravity signal of the seamount is by far larger than the seamount itself, and that the spatial resolution in terms of maximum degree of the gravity field recovered from the satellite mission not necessarily defines the smallest observable object, but it defines the recoverable wavelengths. The integration over the field in the spherical cap shows that a localized mass, which we can call a “bright spot,” does generate a gravity field at a distance of 6° that is still above 10% of its maximum, when considering an expansion up to \(l_{\max }\) = 110 (red curves in Fig. 1). This means that the detection of the field generated by the mass at great distances and at these low degrees depends on the size of the mass and on whether the signal is above the noise level at these low degrees, notwithstanding the diameter of the mass is much smaller than its footprint.

2.3 Illustration of the Spectral Localization with an Example

To illustrate the need and show the correct operation of the spectral localization which supplies a consistent comparison method of the non-stationary signal spectrum and the satellite noise spectrum, we provide an example in Fig. 2: We plot the synthesis of the signal and a noise distribution in the spatial domain, which represent different spatial distributions with identical degree variance spectrum. The function with a regional pattern is the simulated gravity change due to the 2018–2021 volcanic activity at Fani Maoré, the other is a random signal synthesized to have the same spectrum as the former signal, but uniformly distributed over the whole globe. The root-mean-square amplitude of the noise over a patch containing the volcano is much smaller than the volcano signal, showing that if the signal-to-noise ratio amplitude would be a valid criterion, the volcano signal would be well observed, since its amplitude is much greater than the root mean square (RMS) of the noise. The much smaller amplitude of noise is evident in Fig. 2A and B, where the amplitude range of variation of the noise is 0.54 µGal, that of the signal is 8.44 µGal. The difference in amplitude, in contrast with the equal spectrum, is given by the fact that the noise energy is distributed over the entire globe for the global noise, whereas the seamount signal is local, leading to a 16-fold greater peak amplitude of the former.

Fig. 2
figure 2

Example of two global signals with the same spectrum in the SH domain (here expressed as RMS per degree for the first radial derivative of the disturbing potential, \(T_{{\text{r}}}\)). A Seamount gravity signal, using Fani Maoré as an example of a function where most of the energy is concentrated in a small region of the sphere. B One realization of a randomly generated stationary signal, i.e., where variance is homogeneously distributed on the whole sphere. The \(C_{{{\text{l}},{\text{m}}}}\) SH coefficients at each degree \(l\) were randomly sampled, then re-scaled to obtain the same \(\sigma^{2} \left( l \right) = \mathop \sum \limits_{{\text{m}}}^{{}} C_{{{\text{l}},{\text{m}}}}^{2}\) of the SH coefficients of the seamount signal. C The RMS per SH degree spectra of the two SH expansions, which were imposed to be identical. In the map panels A and B, two caps are plotted, referring to the localization windows employed for the localized spectra discussed below: \(\theta_{0} = 6^\circ\) (solid blue line) and \(24^\circ\) (dashed blue line), respectively. A tapered truncation of coefficients centered on SH degree 110 (see Fig. 1) was employed to compute the signals in A and B

Full size image

The spatio-spectral localization is applied on the two example signals, the gravity change due to the seamount and due to the “equivalent noise signal,” which we imposed to be random, stationary and with the same power per SH degree spectrum. The RMS per SH degree spectra obtained by applying spatio-spectral localization on the two example signals are shown in Fig. 3: the gravity change due to a seamount (light blue lines) and due to the “equivalent noise signal” (orange lines). The localization spectrum using the spherical cap radius of 6° (solid line) should be adequate, according to the considerations in chapter 2.2, whereas the radius of 24° (dashed line) unnecessarily dilutes the seamount signal over a greater area, and leads to a downshifted signal spectrum. The noise spectral curve is the same independently of the choice of the cap radius (orange lines in Fig. 3). The signal spectrum with no localization applied is shown with the black dash-dotted line. In contrast, the equality of the global signal and the noise spectra would lead to the false conclusion, that the signal is unobservable, because its spectrum equals that of the noise. The spatio-spectral localization allows to obtain a signal spectrum which can be directly compared to the noise spectrum, and determine the spectral window in which it is above the noise. As further demonstration that the localized spectrum of the signal adequately represents detectability, we calculate the average square root deviation of the signal and noise from zero level, and calculate the ratio between the two values (signal = 3.3 10−11, noise = 1.7 10−12, ratio = 19.4). We can see that the localized spectrum of the signal with the 6° radius results in a spectrum much higher than the noise spectrum, which correctly reflects the almost 20-fold in-cap RMS ratio mentioned above. Here, the given values are a-dimensional, corresponding to the potential divided by the factor GM/r.

Fig. 3 
figure 3

Signal and noise spectral curves: RMS per SH degree spectra obtained by applying spatio-spectral localization on the two example signals shown in Fig. 2: the gravity change due to a seamount (light blue lines) and due to the “equivalent noise signal” (orange lines), which we imposed to be random, stationary and with the same power per SH degree spectrum. The effect of two spherical caps of different radii (6° and 24°) is shown, with solid and dashed lines, respectively—the caps footprint is plotted in Fig. 2A and B. For reference, the spectrum with no localization applied is shown, using a black dash-dotted line. The non-localized spectra of the signal and the equivalent noise are identical, by definition

Full size image

The spectral localization procedure allows the comparison of the spectrum of the local signal, which is essentially zero outside a given spherical cap, with the spectral noise curves of the satellite. As we have illustrated here, the comparison of the localized spectrum of the seamount with the noise spectrum of the satellite is a valid method to determine the detectability. Compared to the signal-to-noise ratio in spatial domain, the spectral approach has the advantage to determine the useful bandwidth of the satellite-recovered signal in terms of detecting the seamount. For a noise curve that increases with degree and order, the useful field should be limited to those degrees in which the noise curve is below the signal spectrum.

3 Review of Seamount Eruptions and Mass Change Estimates for the Selected Seamounts

Here, we provide a review of the volcanic unrests and associated mass and island changes of the volcanoes in our database. The subaerial evolution of the HTHH above sea level is illustrated with the support of satellite Sentinel-1 and Sentinel-2 imagery. Based on the database of mass changes of the volcanoes, we then evaluate the sensitivity of MOCAST+ to detect the volcanic unrests under the most favorable condition that the erupted magma body is replaced inside the volcano by the uprising magma. This condition is acceptable for active volcanoes, in which the magma conduits are pressurized by the uprising magma fed by the crustal or mantle source of molten rocks (e.g., Nicolas 1986; Schmincke 2010). In the less favorable case of a magma reservoir that erupts the magma but does not have a magma replacement and reduces its rock-filled volume, our estimate is an upper limit of the expected mass change. We may argue that in the case of strato volcanoes, the eruptions are successfully increasing the height of the volcano, so that a periodic eruption and recession seems improbable, and the continuous pressurization seems adequate. The image of a magma chamber which remains void before being refilled requires the condition of a rigid surrounding rock able to sustain the local lithostatic pressure and prevent lava and fluid melts to diffuse toward the chamber and refill it. Considering the temperature of the lava above 1000 °C, the temperature alteration affects a region which is considerably greater than the chamber, with partial melting affecting the rocks in decreasing percentage at greater distances. The high temperature of the rocks is incompatible with a rigid surrounding structure, and rather, we must picture the melts to be drawn into the chamber, and the partially molten and viscous rocks will be mobilized to refill the chamber both laterally and from below. Another aspect to be considered is the overpressure at the magma chamber, caused by the reduced density of the hot magma column above the chamber, compared to the thermally undisturbed crustal and lithospheric column: The lithostatic pressure below the hot column is lower by the amount of the density difference integrated along the column and multiplied by the gravity acceleration value. The consequence is that magma is pushed into the magma chamber, leading to a continuous mechanism of its refilling and replenishment of the magma conduits to the volcano summit. According to the diffusion time of the rocks, during the volcanic emission, a mass deficit below the eruption center can appear, but it is recovered driving mass from a much greater area. The overpressure at the magma chamber is increasingly important the deeper the magma source is. The cases we are discussing have magma sources placed at the mid crust or deeper (see details and references below), where the conditions we have just described are valid. This is different for instance to the Kilaulea eruptions, which have a shallow magma chamber of only a couple of km (Poland and Carbone 2018).

We proceed to document the seamount eruptions for the collection of volcanic unrests, and in the end calculate the mass changes for the database of selected seamounts, that are distributed worldwide, as shown in Fig. 4.

Fig. 4
figure 4

Location of the seamounts for which the gravity field is calculated. The volcano Fani Maoré is east of Mayotte. KJ: Kick’Em Jenny Volcano, HTTH: Hunga Tonga Hunga Ha’apai islands and volcano. The other seamounts are labelled with their full name. The background dots (light red) are the seamounts in the Global Seamount Database (Wessel et al. 2010)

Full size image

Unrest in submarine volcanoes may manifest above sea level, leading to the formation of new oceanic islands, a reduction in their surface area to the point of disappearance or changes in bathymetry. This poses a risk due to the emergence of uncharted shallow hazards. We give a detailed description of two big marine volcanoes, of a size class that could be observed through gravity in the future, the HTHH and Fani Maoré volcanoes and further describe in less detail minor volcanoes, limited to submarine extrusions, for which documentation is available.

3.1 Hunga Tonga Hunga Ha’Apai Volcano

The Hunga Tonga Hunga Ha’Apai, or abbreviated Hunga volcano, Kingdom of Tonga, belongs to the Tonga-Kermadec intra-oceanic volcanic arc, generated by the fast subduction of the Pacific Plate. The Hunga volcano emerges 1500 m tall above the sea floor, on a crust 20 km thick (Contreras-Reyes et al. 2011), and topped with an active caldera. The islands Hunga Tonga (HT) and Hunga Ha’Apai (HH) emerge from the caldera rim above the sea level. The stratigraphy of the islands is lava including basaltic andesite to andesite dykes at the base, and above it three sequences of ignimbrites of different ages, topped by volcaniclastic deposits. The historically known eruptions start with 1040–1180 AD, documented in the two HTHH islands. According to Brenna et al. (2022), these sequences are above two or more older volcanic deposits. The successive eruptions of 1912 and 1937 occurred between the two HTHH islands, therefore originating from the caldera rim, from a site at 20°34’S, 175°22’W, 3 km SE from HH. A smaller eruption which did not emerge to the surface was documented for June 1–3, 1988, originating in shallow water 1-km SSE of HH. Lava had erupted from three vents in SW-NE direction extending 100–200 m, but there was no evidence of new island formation in the Annual report of The World Volcanic Eruptions of 1988 (Annual report of The World Volcanic Eruptions 1991). The next eruption on March 17, 2009, emerged from two vents on HH, on the NW coast, and 100-m southwards from the southern tip of the island. The two vents added land made of tephra to the island, increasing the size from 0.51 km2 to 1.42 km2, during the days from March 17 to 21, 2009. Anomalous coloring of the ocean was interpreted as tephra and hydrothermal mineral precipitation in the water column. The tephra is an easily erodible material, which led to a reduction of the island surface to 1 km2 by November 2009. The subaerial erupted volume was estimated to be greater than 0.00176 km3, excluding the submarine deposit, since no information was available on the amount of mass emerging from the caldera to the new island extension. The estimates were based on the remote sensing observations from the instruments ASTER and MODIS, payloads of the Terra and ACQUA satellites (Vaughan and Webley 2010). Next to the new island land, pumice rafts on the ocean were observed also at 40 km from the vent.

The next eruption occurred in 2014 and connected the two islands Hunga Tonga and Hunga Ha’apai through an island bridge. The eruption was active between September 19, 2014, and January 2015 (Garvin et al. 2018). Its evolution has been studied through sea-borne data acquisitions and geologic sampling (Garvin et al. 2018). The volcanic island emerged following the eruption 300 m above the caldera floor. The eruption creating the island has a much smaller dimension than the caldera, occupying only a small portion of the entire caldera.

With the aim to support the understanding of the processes connected to the marine volcanoes, we consider that showing the HTHH evolution in images is very useful. The increase and decrease in the island’s shorelines seen in the images demonstrate that these must be due to vertical movements which remain undetected by the images but in principle are detectable through the associated gravity changes.

3.1.1 Documentation of Eruption History from Satellite Imaging

For the 2014/2015 eruption, the available satellite images document the growth of the island bridge of the Hunga volcano, Kingdom of Tonga, which was destroyed in January 2022, associated to a very strong eruption with a tens of km high ash plume and tsunami (Borrero et al. 2023).

We first describe the observations from the SAR images of Sentinel 1, as shown in Fig. 5, which represents VV polarization in descending mode. We show three steps in the evolution as gray scale images of SAR amplitudes, starting with the situation before the unrest (November 29, 2014), an intermediate stage where the volcano has emerged slightly above the sea surface (December 23, 2014), and after volcanic stability has been reached (February 9, March 5 and 29, 2015). Between HT and HH, the presence of the underlying seamount is not discernable from the surrounding ocean up to the date of November 29, 2014. With the image of December 23, 2014, the appearance of a small cone is visible, with reflectance similar to that of the two islands HT and HH. The HTHH island has grown further on the acquisition of January 16, 2015 (not shown), but shows a strongly reflecting patch on its western side, leading to a homogeneously black patch. The black patch indicates that no portion of the Radar pulse is reflected back to the antenna, which can be due to either a totally absorbing texture, or a well reflecting surface with low backscatter in the direction of the antenna. If the ocean surface is locally covered with a layer of light volcanic material, then the condition of specular reflectivity would be given, or else the patch could indicate the phase of lava outflow. It is not probable it is shadowing, since in the next images, the same flank gives a good backscatter coefficient. The next image of February 9, 2015, shows the fully emerged HTHH island, where though the eastern central flank of the cone has collapsed, and is filled with ocean water. The next image of March 5, 2015, documents a further slight areal expansion of the island, in particular what concerns its connection to the eastern Hunga Tonga island. Starting with the next image of March 29, 2015, the erosional process starts with a visible small reduction of the outline of the island. To define shorelines, the outcome of the edge detector is superposed on the images, and it can be seen that the shoreline of the bridge island is well marked, as is the evolution time. The HT and HH islands show also a few inner features, that are marked as edges of topography or slope change. The image of the first appearance of the bridge island is shown with its outline starting from March 5, 2015. The new island will remain in place until the next unrest in December 2021, demonstrating that the erupted magma was replaced in the caldera by new magma, because the shorelines remained in place, unlikely in the case of a chamber collapse. The chamber collapse would inevitably lead to an island subsidence, which is not compatible with the stable shorelines seen in the satellite images. Only on a much longer time frame of 6 years, its outlines have reduced up to December 2021, after which they increase again. The reduction and later increase in the area above sea level suggests that the island has subsided and then uplifted again, in preparation of the catastrophic explosion of January 15, 2022.

Fig. 5
figure 5

Documentation of the HTHH islands volcanic evolution up to volcanic eruption during the unrest of 2014/2015. SAR Sentinel 1 and optical Sentinel 2 acquisitions. The stability of the outlines of the two peripherical islands suggests that a caldera collapse has not occurred following the eruption that created the bridging island. This observation is important for the mass change estimate

Full size image

Regarding the next unrest of 2021/2022, the images document the preparation phase preceding the January 15, 2022, explosion by 1 month, and the remaining HT and HH islands after the explosion. The subaerial observations of the volcanic unrest leave the submarine changes undiscovered and suggest the usefulness of a combined approach with the gravity monitoring.

The unrest of 2021/2022 starts on December 20, 2021, with a violent eruption, followed by several material emissions in December and first half of January, until on January 15, 2022 (04:14 UTC), a violent explosion occurred generating an eruption column close to 57 km high (Proud et al. 2022) and a tsunami (Crafford and Venzke 2022; Omira et al. 2022; Borrero et al. 2023). It is unknown whether the erupted mass was dispersed as ash in the tens of km high emission, or partly returned as heavier bombs to the volcano flanks. Before the unrest episode starts, the last SAR image is from December 10, 2021, showing the HTHH islands resembling almost the shape they had in 2015, with a minor reduction of the central bridge connecting HT and HH. On December 22, 2021, a new volcano mouth appears NE of the previous cone, with the areal extent of the bridge having almost doubled in surface. Toward the west, the image is disturbed, probably due to the volcanic plume. The next image on January 3, 2022, confirms the enlarged bridge, with the HT and HH islands not having changed. The volcano mouth of 2015 is still present, as is the new vent of December 2021, but with reduced areal extension. It could be that part of the island is submerged, and this is the reason for an altered ocean reflectivity. On the next image of January 15, 2022 (UTC 17:08:36), the central bridge is completely missing, and the two HT and HH islands have reduced size. The next images confirm the stability of the peripheral islands.

For the optical images, we depend on cloud and vaporless days. Sentinel 2 gives a good picture of the bridge and volcanic mouth in 2021 on December 8, 2021. Although very disturbed on the image of December 23, 2021, the bridge is much bigger, and can be seen well on January 2, 2022, with the new volcano mouth north of the previous one on the enlarged bridge. This situation is maintained up to the image on date January 12, 2022. The next image, of January 17, 2022, was taken after the explosion and shows the missing bridge between the reduced HT and HH islands.

Adding the outlines of the islands from the SAR images from 2015 to the optical image of January 2, 2022, suggests that the enlargement of the bridge and to a lesser extent of the HT and HH islands is due to an uplift of the HTHH system. This is suggested by the fact that the 2015 outlines are just a bit smaller for the HT and HH islands, the old volcano mouth (seen well in the image of December 8, 2021) has remained at the same place since 2015, and the southern part of the island appears larger because it has emerged above the height of 2015, putting to light the previously submerged parts of the island.

3.1.2 Mass Change of the HTHH

We estimate first the mass changes for the HTHH unrest of 2014/2015, when the eruption created the bridge between the two HT and HH islands (Fig. 5), leading to an erupted volume of about 90% below sea level emerging above the caldera rim, and 10% forming the subaerial new island (Garvin et al. 2018). It could be argued that the mass estimate is overestimated due to a possible mere shift of the magma from the magma chamber to the surface, assuming the chamber remaining void and without magma replacement from below. This hypothesis would though inevitably lead to a subsidence of the pre-eruption islands, which are placed above the magma source as well. The optical images though document that there is not shore-outline change on the islands before and after the eruption, which indicates that the chamber was refilled from below, remained pressurized, and that therefore the calculation of an extra mass on the top of the seamount, and positioned between the pre-eruption islands is justified.

The volume change due to the eruption of January 2022 (Figs. 6 and 7) has been documented by a post-eruption bathymetry sounding and first estimated to be equal to 6.5 km3 (Stern 2022), with a lowering of the caldera floor to a depth of 800 m from the surface, starting from a depth of 155 m below the sea level. Personal communication by Shane J. Cronin in February 2023 corrected the volume change to 7.9 km3 with an uncertainty of 0.5 km3. There is the uncertainty on the geometry of the mass change, but for our purpose of estimating the gravity change, we find that the exact geometry is not influential, and an extruded cap posed on the summit of the seamount could be used equivalently. For the modeling purpose, we keep to the conservative published value of 6.5 km3 (Stern 2022). To our knowledge, the percentage of mass that was dispersed into the atmosphere and the one that turned back to the volcano close to the volcano is unknown. The volume change of 7.9 km3 though is a net variation, which compares the pre-eruption bathymetry with the post-eruption bathymetry altered by both the explosion and the rocks falling back on to the bathymetry. If we hypothesize the implausible worst case scenario of an explosion with mass redistributed over a much larger radius than the one covered by the ship survey, the spectrum of the gravity field would be reduced at increasingly higher spectral degrees (see Fig. S1).

Fig. 6
figure 6

Evolution of the HTHH islands in the weeks leading to the volcanic eruption during the unrest of 2021/2022. Sentinel 1 SAR acquisitions. The position of the volcano mouth of 2015 remained in place throughout the 2021/2022 unrest, up to the January 15, 2022, explosion that tore away the central island of HTHH

Full size image
Fig. 7
figure 7

Documentation of the HTHH islands volcanic evolution up to the volcanic eruption during the unrest of 2021/2022. Sentinel 2 optical acquisitions in real color (bands red, green, and blue). For comparison in white, the outline of the SAR acquisition of March 29, 2015, is overlaid. The position of the volcano mouth of 2015 remained in place throughout the 2021/2022 unrest, up to the January 15, 2022, explosion that tore away the central island of HTHH

Full size image

3.2 Fani Maoré Submarine Volcano

The Fani Maoré submarine volcano, located 50 km east of the Mayotte volcano, formed during an unrest between 2018 and 2021, and is an example showing that remote gravity observations bring relevant information on eruption monitoring, because the volcanic eruption occurred without evident explosive character, but being accompanied by seismic activity (Berthod et al. 2023; Verdurme et al. 2023). The Mayotte island is located in the Comoros archipelago and is the top of a submarine volcano (latitude −12.8° and longitude 45.15°). A submarine eruption was noticed through onset of a seismic swarm off-centered from the central axis of the Mayotte volcano by 50 km, near to the foot of the volcano, which peaked in June 2018, decaying quickly during 2018, flaring twice in 2019 but with 10 times smaller intensity, reaching low activity from end 2019 onwards. Anomalous ground deformation on the Mayotte island observed by a GNSS network measured eastward horizontal displacements between 21 and 25 cm and subsidence of 10–19 cm. These data together with ocean bottom pressure gauges were inverted to propose magmatic origin of the seismicity and of the ground deformation. The deformation was ascribed to a deflation of a magmatic chamber located at 28-km northwest of the new volcano, at a depth of 24 km (Peltier et al. 2023). Between May 2 and 18, 2019, a multibeam oceanographic campaign discovered a new volcanic cone of 820-m elevation and 5 km3 volume. The acquisition also noted fluid plumes rising from the cone, but not reaching the sea surface. The magma volume erupted from the cone rose from the 3.5-km deep sea floor. Successive campaigns between May 18, 2019, and August 21, 2020, revealed ongoing lava flow at the rate of up to 150–1200 m3/s between July 2018 and June 2019, with an estimated production of another 0.58 km3 of lava up to August 2019. The final estimate of the volume of Fani Maoré is 6.55 km3 (Verdurme et al. 2023). The magma is thought to derive from a primary source at 80–100 km depth which is stored in a magma reservoir in the lithospheric mantle at a depth between 35 and 48 km, and which has a large volume (> = 10 km3) (Berthod et al. 2021).

For Fani Maoré, the composition of the volcanic erupted rocks has been divided into three phases: phase 1 during the 1st year eruption as basanitic magma (4–5 wt% MgO) directly coming from the deep source, phase 2 (starting May 2019) of tephri-phonolitic magma from the base of the crust with more evolved composition (3.4–4.3 wt% MgO) and phase 3, attributed to the exhaustion of the shallow reservoir at the base of the crust, and return to ascent of the magma from the deeper reservoir, returning to composition similar to phase 1 (Berthod et al. 2023). Viscosity was determined on samples, showing that basanite lavas are more prone to silently erupt effusively, whereas phonolite melts can erupt either explosively or effusively (Verdurme et al. 2023). The grain density was calculated for dredged samples from the chemical composition, which leads to an average grain rock density of 2076 kg/m3. Furthermore, density was also measured for the rock samples, with an average sample density of 1900 kg/m3. This density is lower than the grain density due to the porosity; the inferred porosity from the difference of calculated and observed density on average is 30%. The grain density is about 2760 kg/m3, and the dry porous samples have average density of 1900 kg/m3. With the same porosity values, but assuming pores are water-filled, the average density of the wet samples on average is \(\overline{\rho }_{{{\text{wet}}}} \; = \;2760 \;{\text{kg}}\;{\text{m}}^{ - 3}\cdot \left( {1\; - \;0.3} \right)\; + \;1000 \;{\text{kg}}\;{\text{m}}^{ - 3} \cdot 0.3\; = \;2232\; {\text{kg}}\;{\text{m}}^{ - 3}\). With this density value, and a volume of 6.55 km3, the mass change is \(M\; = \;\left( {2232\; - \;1000} \right) \cdot 6.55 \cdot 10^{9} \;{\text{ kg}}\; = \;8.07\; {\text{Gt}}\).

In the case of zero connectivity of the porous voids, with the rock porosity being air-filled, the mass would be 20% lower, equal to \(M_{{{\text{af}}}} \; = \;\left( {1900\; - \;1000} \right) \cdot 6.55 \cdot 10^{9} \;{\text{ kg}}\; = \;5.90\;{\text{ Gt}}\). A global compilation of the relation between connectivity and porosity for effusive volcanic rocks shows that at porosities of 30% connectivity is between 80 and 100% (Colombier et al. 2017). Assuming the value of 80% as the lowest limit, the lower limit for the average expected rock density is \(\rho_{{0.8{\text{wet}}}} \; = \;2760 \;{\text{kg}}\;{\text{m}}^{ - 3}\cdot \left( {1 - 0.3} \right)\; + \;1000 \;{\text{kg}}\;{\text{m}}^{ - 3}\cdot 0.3 \cdot 0.8\; = \;2172\;{\text{ kg }}\;{\text{m}}^{ - 3}\). The emplaced mass change for the 80% connectivity volcanic rocks, with partially air filled porosity, and partially water-filled porosity, would be \(M_{{0.8 {\text{wf}}}} \; = \;\left( {2172\; - \;1000} \right) \cdot 6.55 \cdot 10^{9} \;{\text{kg}}\; = \;7.68\;{\text{ Gt}}\). It can be assumed that the porosity will be reduced toward the central axis of the volcano, depending on the thickness of the overburden, so this estimate can be considered a conservative value.

The replacement of the erupted magma in the magma chamber from below is not known, so the net mass change in theory could be smaller. Considering the magma chamber being at a depth of 35 km, it seems mechanically impossible that a hollow chamber remains, but rather that it is at least partially refilled with magma due to the high surrounding lithostatic pressure. The less the chamber is refilled, the smaller will the net mass change be.

3.3 Surtsey Island

Surtsey island is located SW of Iceland and belongs to a volcanic system that probably started activity 100,000 years ago. The volcanic system is separated from Iceland and includes 17 volcanoes emerging from the sea floor to the surface and another four submarine volcanoes. Further hills and peaks on the seafloor are thought to be remnants of more than 40 submarine late Pleistocene volcanoes. In historic times, Surtsey was active in 1963–1967, eruption which we include in our review.

The unrest of Surtsey during 1963–1967 lead to the formation of the Surtsey island, rising from the sea floor, with an erupted volume of 1.1 km3, consisting of 70% tephra and 30% basaltic lava. The subaerial tephra is estimated to have 45–50% porosity. The volcano measures 13.6 km2 at its base, narrowing to the area of Surtsey island presently of 1.41 km2. The volcano is maximum 2.9 km wide. The height of the volcano is 285 m, measured from the sea floor 130 m below sea level. The tephra was formed from basaltic volcanic ash, produced by quenching of hot magma by cold sea water. Surtsey has given the name to this rock type, described here for the first time, and giving it the name “Surtseyan tephra” (Walker and Croasdale 1971). Following the explosive tephra phase of the eruption, the effusive phase deposited successive basalt shields between 1964 and 1967, with volume of 0.4 km3, reaching 170-m thickness above sea level. A morphological feature that is interesting in the context of the area increase in HTHH in December 2021 is that the lava forms submarine deltas in Surtsey, in front of the shore, with shallow slopes. This delta is made of brecciated lava, which reaches 130-m thickness.

To model the geometry of Surtsey, we approximated it by a cone of 13.6 km2 at its base, 1.4 km2 at sea level and 155 m tall from the bottom to the sea surface. The submarine volume is VSM = 1.16 km3 for a cone of these values. The equivalent area of the submarine part, calculated assuming a constant area, is equal to \(A_{{{\text{SAEQ}}}} \; = \;\frac{1.16}{{0.155}}\;{\text{km}}^{2} \; = \;7.48\;{\text{ km}}^{2}\). The subaerial volume, in the approximation of a cylinder, is \(V_{{{\text{SA}}}} \; = \;1.4 \cdot 10^{6} \;{\text{ m}}^{2} \cdot 170\;{\text{ m}}\; = \;0.24\;{\text{ km}}^{3}\). The ratio of the subaerial part to the total volume is \(r_{{{\text{SA}}}} \; = \;\frac{{V_{{{\text{SA}}}} }}{{\left( {V_{{{\text{SA}}}} \; + \;V_{{{\text{SM}}}} } \right)}}\; = \;0.17\). The subaereal density of tephra is \(\rho_{{{\text{tephra}}}} \; = \;2800 {\text{ kg}}\;{\text{m}}^{ - 3} \cdot \left( {1\; - \;0.5} \right)\; = \;1400\;{\text{ kg}}\;{\text{m}}^{ - 3}\), assuming 50% air-filled porosity. The submarine density with water-filled porosity is \(\rho_{{{\text{tephra}}}} \; = \;2800 {\text{ kg }}\;{\text{m}}^{ - 3} \cdot \left( {1\; - \;0.5} \right)\; + \;1000 {\text{ kg }}\;{\text{m}}^{ - 3} \cdot 0.5\; = \;1900\;{\text{ kg }}\;{\text{m}}^{ - 3}\). We now calculate the rock density for a given volume of a given percentage of basalt and tephra: subaerial rock density of 30% basalt and 70% tephra \(\rho_{{{\text{rock}},\;{\text{SA}}}} \; = \;2800 {\text{ kg }}\;{\text{m}}^{ - 3} \cdot 0.3\; + \;1400{\text{ kg }}\;{\text{m}}^{ - 3} \cdot 0.7\; = \;1820\;{\text{ kg }}\;{\text{m}}^{ - 3}\); the same rock in submarine position and full connectivity of pores, and considering the density change against water is \(\rho_{{{\text{rock}},\;{\text{SM}}}} \; = \; {2800 {\text{ kg }}\;{\text{m}}^{ - 3} \cdot 0.3\; + \;1900 \;{\text{ kg}}\;{\text{ m}}^{ - 3} \; \cdot 0.7}- \;1000\;{\text{ kg}}\;{\text{ m}}^{ - 3} \; = \;1170\;{\text{ kg}}\;{\text{ m}}^{ - 3}\). With the volumes calculated before, the subaerial mass change is \(M_{{{\text{SA}}}} \; = \;0.24\;{\text{ km}}^{3} \cdot 1820\;{\text{ kg}}\;{\text{ m}}^{ - 3} \; = \;0.437\;{\text{ Gt}}\), and the submarine mass change is \(M_{{{\text{SM}}}} \; = \;1.16\;{\text{ km}}^{3} \cdot 1170\;{\text{ kg}}\;{\text{ m}}^{ - 3} \; = \;1.357\;{\text{ Gt}}\).

3.4 Minor Submarine Volcanic Eruptions

The minor submarine volcanic eruptions refer to the Axial, El Hierro, Loihi and Kick 'em Jenny seamounts. The Axial seamount (USA) is the most active underwater volcano in the NE Pacific, located in the Juan de Fuca Ridge and above the Cobb hot spot (Chadwick et al. 2016). More than 50 eruptions have occurred in the past 1600 years (Clague et al. 2013), the latest events on record in 1998, 2011 and 2015. The 2015 eruption is of note, being the first one observed in real-time by the Ocean Observatories Initiative (OOI) instrument network. Time-lapse bathymetry between 2015 (after the eruption) and two different prior surveys (2006, 2013) allowed to identify the area interested by lava flows: three northern flows with maximum thickness ranging from 67 to 128 m and eight smaller southern flows 5–19 m thick (Kelley et al. 2014; Chadwick et al. 2016). The northern flow represents ca. 92% of the total erupted volume, which was estimated to \(1.48 \cdot 10^{8} {\text{m}}^{3}\). As a term of comparison, the estimate for the 1998 and 2011 eruptions amounts to \(0.31\) and \(0.99 \cdot 10^{8} {\text{m}}^{3}\), respectively.

The correlation between the geochemistry of the 2015 lava flows and their location, together with the rate of inferred inflation of the caldera since 2011 (Nooner and Chadwick 2016), suggests a recent magmatic recharge and a decrease in magma residence time. Chadwick et al. (2016) observed that Axial could have been evolving toward a period of larger and more frequent eruptions. However, the monitoring made since 2015 has shown that magma resupply, which allowed Axial to re-inflate to 85–90% of the deflation it experienced in 2015, has since slowed down, consistent with already observed decadal variations in the magma supply rate (Chadwick et al. 2022).

The El Hierro submarine eruption occurred at the foot of the Canary island El Hierro, between October 2011 and February 2012. The eruption was non-explosive and in successive lava flows created a seamount 325 m tall, reaching 88 m below sea level (Rivera et al. 2013). The bathymetry change is well documented through monthly multibeam ship sounding sent to monitor the volcano during evacuation of parts of the El Hierro island (Carracedo et al. 2015). A strong increase in seismicity anticipated the eruption starting from July 2011, occupying a seismic volume (10 km by 30 km horizontally, by 15 km in depth) by far larger than the seamount diameter (1 km). The seismicity migrated from North to South of the island over a distance of 30 km, at depths between 10 and 25 km, in the mantle (Carracedo et al. 2015). One GPS station on the island recorded inflation before the eruption, followed by an episode of deflation during the early stages of the eruption, followed by a return to inflation with net and stable uplift of 22 cm after the eruption and unrests ended. The repeated bathymetry measurements demonstrated a net increase in mass accumulated at the volcano, with a stable bathymetry at tens of kilometers from the volcano, which does not fit the model of a caldera emptying and collapsing due to release of magma. This observation suggests that our mass change is correct, since the magma beneath the volcano is drawn from a large and deep area, hardly influencing the local gravity field. El Hierro is another example of a seamount growth exceeding 300 m, which would have remained unrecognized, if it had not been specifically visited by the scientific vessel sparked by the appearance of anomalous seismicity and the closeness to the inhabited island. The volume is documented, the rock is basalt, but the porosity is not known. Therefore, the mass estimate is an upper limit, which could be lowered by up to a factor 30% in case of levels of porosity as those dredged at Fani Maoré.

The Loihi Seamount (USA) is located near the island of Hawaii and due to its shallow summit depth of about 1 km could evolve in the future into a new island of the archipelago (Lipman and Calvert 2013). The Kick 'em Jenny Seamount (Grenada) is an underwater volcano positioned under an important navigation route, as it is close to the island of Grenada in the Caribbean (Carey et al. 2016); there have been at least 12 eruptions since its discovery in 1939. However, the volume erupted by Kick 'em Jenny is very small, so it is expected to be below the threshold of what can be observed by the gravimetric satellite. The documentation of erupted volumes are lacking, so our mass value is an example estimate and considerably greater volumes could occur. 

For these minor submarine volcanic eruptions, we refer to the publications mentioned above, in which the erupted volumes are published. The eruptions are submarine, so the density contrast must be calculated against water, and we assume the volcanic eruption consisting of basalt of density 2800 kg/m3, thus a density contrast of 1800 kg/m3.

4 Summary of the Mass Changes of Submarine Volcanoes

In Table 1, the volumes of the mass change due to an eruption, the area of the volcano summit and the corresponding thickness of the eruption placed at the summit are given for the review of selected seamounts. We recall that the mass change is distributed on the summit. The density of the mass surplus is for a basalt eruption, taken as standard except for the HTHH, Fani Maoré and Surtsey, where a lower density of the erupted rocks has partly been documented. The mass changes based on the bathymetry changes at the volcano are an upper limit of the net local mass change deduced from the bathymetry change, because the replacement mechanism at the magma reservoir is not known and not defined.

Table 1 Parameters used to model the synthetic eruptions for the seven seamounts
Full size table

For each of the volcanoes, we calculate the gravity field, with the expression of an extruded spherical cap set at the top of the seamount. Next to the geographical gravity field, we calculate the SH expansion of the field. As explained before, the spectral localization is needed to obtain the localized spectral signal curves to be compared with the noise curves.

5 Sensitivity Evaluation of MOCAST+ to Seamount Eruptions

Having calculated the gravity signals and the SH expansion of the signals of our seamounts database, we compare them with the mission error spectral curves. We consider the error degree curves of different MOCAST+ configurations, given in Table 2, which serves as a guide for the legend entries of all the spectra shown hereinafter (Migliaccio et al. 2023). The payload consists in a quantum gradiometer and clock, on a satellite which flies in formation on polar and inclined orbits, called Bender formation (Bender et al. 2008). The gradiometer measures the second derivative of the gravity potential field along the radial (Tzz) direction, or in the direction orthogonal to the orbit plane (Tyy). The gradiometer is a cold atom interferometer with Sr atoms, with a noise level 1 mE Hz−1/2 in the frequency band from 5 × 10–6 to 1 × 10−2 Hz, with steeply increasing noise at higher frequencies. In the original paper also, the gradiometer with Rb atoms is discussed, with slightly higher noise level (Rossi et al. 2023). The attitude compensation is considered to be either at full level (ideal case), or with an accuracy of 1 nrad s−1. The compensative attitude control is necessary due to the rotation of the satellite around the y-axis along the orbit which introduces Coriolis forces on the Tzz component which largely increase the noise. The gravitational potential difference at satellites carrying atomic clocks linked with a laser beam is recovered assuming an optimistic scenario with an error of 0.21 m2 s−2 Hz−1/2. The pairs or triplets of satellites are placed in a polar orbit (inclination 89°) and an inclined orbit (inclination 66°), at mean altitudes of 371 km and 347 km, respectively. The distance between the in-line satellites is 1000 km, which means that for three inline satellites, the leading satellite has 2000-km distance from the third satellite. The EGM2008 gravity potential model was used as static reference model, to which the non-tidal mass variations in atmosphere, ocean, hydrology and solid Earth were added using the ESA Earth system model (Dobslaw et al. 2015, 2016; Rossi et al. 2023), used also for dealiasing purposes. Combinations of payloads and satellites are listed in Table 2, with nomenclature in the configuration description as follows: D is the distance between two satellites, Tzz or Tyy is the gradiometer component measured and optimal noise PSD (power spectral density) for the gradiometers is intended to mean Sr atoms with ideal compensation for drag and rotation. 1 or 0.1 Hz clock indicates the clock observation sampling rate. The best performance is found for the two triplets of Tzz components and clocks, assuming full compensation of rotation. Perfect compensation can be considered as a theoretical limit, with a more realistic case assuming the presence of attitude control error. In that case, the best performance is obtained for five satellites with Tyy, three on inclined, two in polar orbit, one satellite with Tzz on polar orbit and compensation at 1 nrad/s level and clocks on all six satellites. This solution is slightly better than the one with two triplets of Tyy components and clocks, the one included in Table 2. For details on the assumptions on payload and orbits and gravity field simulation results, please refer to Migliaccio et al. (2023) and Rossi et al. (2023).

Table 2 Details on the configurations considered for MOCAST+ and corresponding names in the legend of the following figures
Full size table

All the error estimates were formulated for a timespan of 1 month. In the case of the “MOCAST+ 2 couples nominal” configuration, we were also provided with the results of simulations covering 4 months. Apart from the specific case on short-term detectability, a timespan scaling was applied to account for the reduction in error due to longer observation intervals, as for 1 or more years. This reduction, arising from the improved coverage and stacked measurements, is modeled with the following scaling relationship:

$$\sigma^{2} \left( l \right)_{{\Delta t_{1} }} = \sigma^{2} \left( l \right)_{{\Delta t_{0} }} \cdot \left( {\frac{{\Delta t_{0} }}{{\Delta t_{1} }}} \right)$$
(2)

where \(\sigma^{2} \left( l \right)\) is the error degree variances at the SH degree \(l\) formulated for the solution timespans \(\Delta t_{0}\) (original time interval) and \(\Delta t_{1}\) (time interval we are scaling to).

For comparison, we consider the error degree variances of GRACE (Mayer-Gürr et al. 2018; Kvas et al. 2019). The GRACE error curves employed for the analysis are the average error degree variances of all the 162 monthly solutions of the ITSG-GRACE2018 model (Kvas et al. 2019), resulting in an average error of a monthly solution. It should be noted that this is not equivalent to stacking 162 1-month solutions, which would incorrectly result in a significantly lower error estimate. The use of noise analysis in the SH domain alone is dictated by the fact that the MOCAST+ mission concept is very different from the GRACE double pair concept: For instance, no inter-satellite distance measurement is carried out. Hence, the only coherent comparison method is through the comparison between signal and error degree variances. This does not imply that an alternative local analysis of the Level 0 or Level 1 GRACE data analyzed over the volcano region could not give a lower noise, but this more specialized analysis of the GRACE data is out of the scope of this paper. For instance, the range rate changes or the acceleration changes between the GRACE satellite couple could be analyzed, or an ad hoc regional solution could be obtained trough a parametrization tailored to seamount retrieval, which could be the subject of further work. Our analysis aims at reproducing a seamount detection strategy based on analysis of readily available Level 2 products, such as the monthly gravity field SH expansions provided by different data processing centers to the ICGEM data services (Ince et al. 2019). This setup, albeit partial by design, ensured full consistency with the products of MOCAST+ simulations.

The localized signal spectra compared to the MOCAST+ spectral error curves are shown in Fig. 8 and are seen to be quite low: Five out of seven of the analyzed seamounts are non-detectable in the full degree range. The localized spectrum is obtained for a 6° degree radius cap, the optimal radius according to the trade-off shown in Fig. 1. A larger radius allows to obtain lower degrees in the calculated spectrum, at the cost of reduced spectral amplitudes of the volcanic fields (an example of the effect of using a larger radius is shown in the supplementary material, Fig. S2).

Fig. 8
figure 8

Comparison of the spectra of the submarine volcanic eruptions with respect to the spectral error curves of MOCAST+ , GRACE. A The signal curves of all seamounts compared to error curves of GRACE and the most favorbale MOCAST+ configuration (MOCAST+ 2 triplets Tzz + C in Table 2). B The signal spectral curves of HTHH 2022 and Fani Maorè compared to different configurations of MOCAST+. The signals for HTHH 2022 and Fani Maorè are too small with GRACE, but could be seen by MOCAST+ 

Full size image

For the minor seamounts, the calculated gravity signals of the seamounts are small, so it is a challenge to observe them from satellite with a short time of observation, as 1 month. Increasing the time window of data acquisition to several years, the noise curves are lower, and we find that the bigger eruptions as HTHH of 2022, of Fani Maoré and Surtsey have a chance to be seen by the MOCAST+ mission, still being at the limit to be observed by GRACE.

Considering the different configurations proposed for MOCAST+ mentioned in Table 1 (Rossi et al. 2023; Migliaccio et al. 2023), the configuration with the gradiometer component Tyy is not favorable in detecting the volcano signal, the Tzz component with rotation compensation is required to achieve detectability. Further, a mission relying on only the clock observations would not detect any of the volcanoes.

6 Discussion

This work aims to detect submarine volcanic eruptions with quantum technology instrumentation on board a future satellite mission constellation. To prepare for such endeavor, it is necessary to (1) define the expected mass changes of submarine eruptions, (2) calculate the expected gravity field and (3) compare this field to the noise characteristics of the quantum gravity mission. We have reviewed the documented submarine volcanic eruptions and defined the mass changes for seven cases, which include silent eruptions with no external signature as the Fani Maoré volcano offshore Mayotte, and the massive HTHH eruptions of Hunga at which repeatedly islands are created and destroyed, or silently change their shape due to subsidence or uplift, attributable to the caldera depressurization or inflation. The silent eruptions are very difficult to cover globally through alternative data acquisitions, because the volcanism is distributed over the entire oceans, including extremely remote areas, and because it is difficult to distinguish a swarm seismicity accompanying a magma ascent at an uncharted seamount from regular tectonic seismicity. If the volcano reaches the surface and creates an island, remote sensing with SAR and optical imaging can catch the newborn eruption, but this cannot be done if the seamount does not perch the ocean surface. An example is the Fani Maoré volcano which rose 800 m tall, remaining below the ocean surface, and was only detected because the GPS stations on the Mayotte island recorded a deformation due to a source offshore the main island and an anomalous seismicity. This means that such a silent eruption can be ongoing in remote areas along the mid ocean ridge for instance, without being detected, but potentially creating hazard to ships and submarines, or also to submarine communication cables. It seems, therefore, important to investigate the possibility to detect these mass changes through the gravity field products of the next-generation gravity satellites. This investigation is timely due to the planned NGGM/MAGIC double GRACE-type mission constellation expected to be active by end of the 2020-ies, the quantum technology Pathfinder mission expected by end of the 2020-ies and the quantum gravity technology gravity mission following the Pathfinder explorative mission. The definition of the setup of a satellite mission constellation by a Space Agency requires the Mission Requirements Document driven by the scientific community user needs. We have seen above that the noise curves of the retrieved gravity field depend on the mission constellation, in case of the MOCAST+ mission on the instrument noise level, on the number of satellites in the constellations (2 couples in Bender formation or two triplets in Bender formation), on the distance between the satellites flying in formation (Rossi et al. 2023; Migliaccio et al. 2023). The technological choices are a trade-off between budgetary and scientific optimizations, and one criterion is that the mission must fulfill the threshold observation requirements.

In the study, we have estimated the upper limit of mass changes of documented seamount eruptions including the HTHH Hunga eruption and Fani Maoré, Surtsey and other minor documented eruptions. We find that the large scale events can be observed with the quantum technology gradient satellite mission. The gravity potential field observation is the only means to detect a seamount growth in case the island remains under water and cannot be detected by optical passive or active SAR observations. The mass change induces a change in the potential gravity field and the geoid, which in principle affects the sea level, opening both satellite altimetry and satellite gravity as space borne methods to detect the change.

The comparison of the volcano gravity signal with the noise of the recovered gravity signal of MOCAST+ is done in the spectral domain, because the noise is published in terms of uncertainty on the spherical harmonic expansion coefficients by the partners of the project in charge of calculating the gravity retrieval performance of the MOCAST+ mission (Rossi et al. 2023; Migliaccio et al. 2023). We have explained above that due to the localized nature of the volcano signal, the spectral localization of the volcano SH representation is necessary, to make the two quantities comparable. Given the colored nature of the degree variances of the two spectra, with the signal decaying and the noise increasing with higher SH degrees, there is typically a spectral window for which the volcano signal is above the noise curve. This window defines the segment of the spectrum in which the volcano signal will be detectable by the MOCAST+ mission. The comparison with the noise curve of GRACE-FO demonstrates the improvement that the MOCAST+ mission would allow, which is in particular defined in terms of:

  • Ability to observe a larger number of phenomena (i.e., sensitivity to smaller volcanic volume increases the number of observable eruptions)

  • Increase in the sensitivity to lower mass changes (i.e., volcanic eruptions with smaller erupted volume)

Our main findings are described in Table 3.

Table 3 Overview of the detectability of geophysical signals by MOCAST+ and GRACE considering a standard L2 data product published in terms of Stokes coefficients of the recovered gravity field. With \(T_{{\text{r}}} {\text{Max}}\), we report the absolute peak in the radial derivative of the disturbing potential \(\left( {T_{{\text{r}}} } \right)\) at ground level, when synthesized up to d/o 120
Full size table

The detectability considerations discussed thus far focus on the sensitivity of the mission proposal. Another critical factor is the mitigation of disturbances from ocean currents through models of ocean circulation. If not accurately modeled, these effects may diminish the detectability of the seamount, as the signal needs to be distinguished from the oceanic signal. The ocean contributes to gravity noise, and its correct modeling or averaging is crucial for enhancing detectability. In our work, we assess detectability based on the noise curves of the gravity field retrieval over a specific time interval. The noise estimate is derived from the difference between the input and recovered fields. The input field contains a time-variable component using the atmosphere–ocean time-variable signal, while the recovered field is obtained from satellite observations over a given time interval. This implies that the noise estimate includes errors in retrieving the time-variable gravity signal but not errors in modeling ocean variation. The development of tidal oceanic models relies on hydrodynamic modeling, calculating oceanic phase and amplitude at tidal frequencies, empirical models and tuning hydrodynamic models with observed terrestrial and altimetric data. Examples include models based on satellite altimetry like DTU10 (Cheng and Andersen 2011), the Goddard Space Flight center models (GOT) series (Ray 2013; Ray and Erofeeva 2014), the Ohio States University, OSU models (Fok 2012), the EOT models (Hart-Davis et al. 2021) from the German Geodetic Institute and the FES2014 (Lyard et al. 2021) tidal model which assimilates altimetric time series for a numerical geophysical modeling of the ocean tides. The TPXO series from Oregon State University (Egbert and Erofeeva 2002) provides a user friendly calculation service of ocean tides, calculated models and software. These models are constantly improving, with errors as low as 5–10 cm for remote areas like polar regions (Sun et al. 2022). Modeling non-tidal ocean variations is more challenging, and insights on handling this noise interfering with the detection of tectonic or volcanic transients come from studies of ocean bottom pressure gauges. The recent study of Otsuka et al. (2023) in the Nankai Trough offshore Japan applied Principal Component Analysis to 3 years of ocean bottom pressure gauge data, successfully separating transient tectonic signals from total variability. The challenges faced in analyzing gravity changes in space and time for detecting volcanic mass change signals align closely with those addressed using PCA analysis based on results from ocean bottom pressure gauges. The spectral curves we have obtained for the different volcano eruptions allow to give a feedback to the choices on payload and satellite mission concept of MOCAST+ , which both affect the steepness of the error curves in the spectrum. The improvement in the gravity field estimation at low degrees through the two triplets is not relevant for seamount detection. The two couples configuration is equivalent at the degrees of interest for the detection of the seamounts. The End-to-End (E2E) simulations of gravity missions are made using a synthetic model of the gravity field variation time (AOHIS model) which considers the mass variability in the atmosphere (A), oceans (O), terrestrial water storage (H), continental ice-sheets (I) and the solid Earth component, with the ongoing glacial isostatic adjustment and the Sumatra–Andaman M9.1 Earthquake of 2004 (S) (Dobslaw et al. 2015). Our simulation could render the AOHIS model more complete, with the addition of the seamounts. The lower noise in the gravity field estimation at low degrees results in the possibility to improve the detection of the geophysical phenomena, which contain also low energy at harmonic degree < 25. In any case long period components of hydrology, large-scale atmospheric circulation, ocean tides would definitely benefit from lower errors at these low spatial frequencies.

The MOCAST+ constellation with clocks and gradiometers could lead to the detection of seamount eruptions with volumes in the order of 1.5 km3 ± 0.5 km3 or greater, assuming that the magma chamber be refilled from the mantle source. We mentioned that the calculated mass change is an upper limit, as it does not consider the magma replacement mechanism at the magma chamber. The detection of the evolution in time of the mass change gives a constraint on formulating models for the magma flow into the chamber. The detection of a zero gravity change, for a mass change above the detection limit, would give a limit on a mere mass displacement from the source to the volcano summit, and the probable following mass increase an indication on the chamber refilling with an information on its timing. The Axial seamount is the best monitored seamount presently, where a co-eruption small subsidence of the caldera (2 m), followed by an uplift which started immediately after the eruption, is an indication that the refilling of the chamber started immediately after the eruption (Chadwick et al. 2016 and Axial Seamount Cruise Reports).

7 Conclusions

Our efforts to define a review of documented seamount eruptions, including the mass change, the gravity field change and comparison of the signal spectrum to the noise spectrum of the MOCAST+ mission proposal, were successful in showing that MOCAST+ would be able to detect seamount growth that has occurred silently and unobserved by remote sensing techniques. In this application, MOCAST+ would be an improvement with respect to GRACE-FO, in the aspects that can be summarized as follows.

The gravity monitoring is demonstrated to be able to detect submarine eruptions starting from net mass changes of 7.7 Gt, with a latency of 1 year.

The required latency decreases for larger events, because the signal is above the noise curve after a shorter time of observation. The retrieval error of the gravity field estimates decreases with increasing integration in time and averaging of the observed gravity signals. The volcano eruption being a permanent signal and not a transient, it is valid to integrate the observation over a given time interval, with the consequent noise reduction in the observed fields.

Few alternatives exist to detect silent submarine eruptions—the Fani Maoré is a fitting example: anomalous seismicity and deformation (on the nearby Mayotte volcano, by GNSS) have been detected there. However, only after a ship expedition had been sent, the enormous new 800 m tall seamount was detected. The same eruption far from inhabited islands would not have been identified as such, producing a hazardous uncharted seamount, dangerous for submarines and ships, in case it grew closer to the sea surface.

In addition to our review and detectability assessment, the gravity fields we calculate for the submarine volcanoes are suitable to integrate the time variable gravity products (including atmosphere, ocean, hydrology, isostasy and only one earthquake signal) used in mission simulations and for purposes such as correcting the time-varying gravity field in climate change research using gravity field observations.