Candidates for volcanoes under the ice of Antarctica detected by their gravito-topographic signal

We analysed the newest gravity field models (given in the form of spherical harmonic expansion to degree and order 2190) EIGEN 6C4 and RET 14, which consist basically from the gravity data from EIGEN 6C4 plus the bedrock topography from the Bedmap 2 model, describing the topography of the ground under the ice of Antarctica. From these models we computed the gravity anomalies, the Marussi tensor of the second derivatives of the disturbing potential, the gravity invariants and their specific ratio, the strike angle and the virtual deformations. We applied these results for a selected part of Antarctica; not too far from the Lake Vostok we discovered at least two objects that might be subglacial volcanoes. We present all predictors and arguments we have available to support this hypothesis, but we are well aware that the “last word” waits for other specialists. We provide geographic coordinates on the ice surface where to dig for the


Introduction
We use the best now available gravity data (the harmonic geopotential coefficients from the gravity field model RET 14 [Hirt et al. 2016], focused on Antarctica, transformed to the various gravity functions) and the bedrock topography data [Bedmap 2, Fretwell et al. 2013] to detect possible subglacial volcanoes and other features (like subglacial lakes and rivers). This paper is a continuation of Klokočník et al. [2016]; it deals with possible subglacial volcanoes (glaciovolcanoes) discovered near the Lake Vostok (LV).
Volcanoes are common on this planet and there is about 30 volcanoes (not subglacial) already known in Antarctica, see, e.g., AntarcticGlaciers.org, last updat-ed in March 2014; LeMasurier &Thomson 1990 andLough et al. 2013 for discovery of two active subglacial volcanoes in Marie Byrd Land from seismic observations. Many more such volcanoes can still be hidden under the ice [e.g., Walton 2013, p. 46].
We can only speculate about the number and age of the subglacial volcanoes in Antarctica. Knowledge of their location is important for further investigations on their role in glacial motion in Antarctica as well as seismic activities and the tectonic setting. Their possible role in glaciation and deglaciation of Antarctica is also a highly interesting topic. So we have good reasons to seek for candidates of subglacial volcanoes. Active subglacial volcanoes are known to be located in Iceland, British Columbia and Yukon. What we could expect in Antarctica due to the local conditions (location under the ice) are small slopes of the volcanoes normally associated with lava-dominated shield volcanoes of Hawaiian type.
As for detection of subglacial volcanoes in Antarctica we do not consider seismic tools or magnetic data, but the density contrasts due to a volcano, inducing their specific "gravity signal"; we also use the bedrock topography. Both data types are used simultaneously. We are very well aware of non-uniqueness of the gravity data (one and the same gravity signal observed on surface may mean more than one type of features hidden under the surface and generating the causative density anomaly). Thus, we prepared a list of the gravity symptoms for a huge single mountain or volcano (Section 4).
At LV, a tectonically controlled subglacial lake, the tectonic processes provided the space for a unique habitat and recent minor tectonic activity could have the potential to introduce sufficient amounts of thermal energy into the lake [e.g., Studinger et al. 2003a]; a seismic activity was recorded nearby [Studinger et al. 2003b[Studinger et al. , 2004, thus we speculate that volcanoes might be expected not far away. Our study area (Figure 1a) is around LV and in the west direction to GSM (Gamburtsev Subglacial Mountains).
The core of our method [Klokočník et al. 2016[Klokočník et al. , 2017b is in the use of various gravitational (or gravity) aspects, namely the gravity anomalies or disturbances Δg, the components of the Marussi tensor Γ of the second derivatives T ij of the disturbing potential (especially the radial component T zz ), the gravity invariants I 1 and I 1 , their specific ratio I, the strike angle θ and the virtual deformations vd. Each such gravity aspect contains different information and has a different sensitivity to the underlying causative density variations [see Figure 2 in Klokočník et al. 2016]; the traditional gravity anomalies are already not sufficient, more gravity aspects are used. The theory with many examples (Figures) can be found in Klokočník et al. [2014Klokočník et al. [ , 2016Klokočník et al. [ , 2017b.
The theory underlying our methodology is described mainly in Pedersen and Rasmussen [1990] and Beiki and Pedersen [2010]. Many examples for local use of the gravity aspects can be found in the literature dealing with petroleum, metal, diamond, economic minerals, oil and gas deposits, fault location, groundwater explorations etc. [see, e.g., Saad 2006, Murphy and Dickinson 2009or Mataragio and Kieley 2009. The Marussi tensor provides more complex information than the ordinary gravity anomalies only; T zz informs about the setting (target location), the other components T ij of Γ refer to the orientation of the causative structure. The invariants can be looked upon as non-linear filters enhancing the sources having big volumes; they discriminate major density anomalies into separate units. The ratio I can indicate two dimensionality of the causative body. The condition I=0 is necessary but not a sufficient condition for two-dimensionality. The strike angle θ tells us how the gradiometer measurements (used to derive Γ) rotate within the main directions of the underground structure; when I=0, the values of θ may indicate a dominant 2D structure [used to predict possible oilfields in Klokočník and Kostelecký 2015]. The virtual deformation vd is an analogy with the tidal deformation and characterizes the "tensions" (compression and dilatation) generated by the causative body; one can imagine the directions of such a deformation due to "erosion" brought about solely by the "gravity origin" Klokočník et al. 2013, p. 90]. The vd provide a dynamical information although they are, as well as all the gravity aspects, computed from the static gravitational geopotential.
The theory used (with our own contribution to it), the method and the notation (with various results for different parts of the world) can be found in Kalvoda et al. [2013], Klokočník and Kostelecký [2015] or Klokočník et al. [2010Klokočník et al. [ , 2014Klokočník et al. [ , 2016Klokočník et al. [ , 2017a. The computation of all the aspects listed above is not trivial. It was done by means of software developed mainly by B. Bucha [see Bucha and Janák 2013, where reader can find technical details] and by our own software.
The input to the computations is always a set of harmonic geopotential coefficients of the given static gravity field model (here RET 14), the outputs are the individual gravity aspects, which are then plotted (individually with non-linear scales) with a 5x5 arcmin step in latitude and longitude. We do not use any filter during our computations. About plotting and coordinate systems used see Klokočník et al. [2017b], Sections 2.2 and 3.5.
We have scanned the whole world by means of several recent gravity field models with a step of 5x5 arcmin in geodetic latitude and longitude, including EIGEN 6C4 [Foerste et al. 2014]. This 5 arcmin step corresponds to the resolution of EIGEN 6C4 (and RET 14, too, see Section 2.1) which is dictated by the spherical harmonic expansion of the gravity model available to degree and order 2190. We always use this full resolution. The RET 14 is intended only for Antarctica, where it is valid (Section 2.1); for many examples see Kalvoda et al. 2013, Klokočník et al. 2014, 2017. We work also with the bedrock topography; it is given in a network 1x1 km (in a specific x,y coordinate frame of Bedmap 2), but it does not mean that everywhere the topography data have actually such a resolution (more about this topic is in Sections 2.1 and 2.2).

The data files
In general terms, the "gravity" is represented by a static model of gravity/gravitational field of the Earth in terms of the harmonic geopotential coefficients (also known as Stokes parameters) to a certain degree l and order m, C lm , S lm , derived from a combination of satellite and terrestrial data. The "topography" is represented by a model of the bedrock topography for Antarctica achieved dominantly via the ice penetrating airborne radar (also called RES, radio echo sounding) data.
In particular, we use the best data now available for Antarctica, i.e. RET 14 [Hirt et al. 2016]. It is a degree 2190 gravity field model SatGravRET2014 (shortly RET 14), given as a set of harmonic geopotential coefficients, meaningful however only for the continent of Antarctica (not for the whole world). Roughly speaking it combines the global gravity field model EIGEN 6C4 and the Bedmap 2 topography (see below). More precisely: the RET 14 combines the predecessors of EIGEN 6C4 with other and very important gravity data sets; it combines the ITG-GRACE2010s and the unconstrained GOCE TIM5 satellite gravity models to degree and order 180. These models describe the long-and medium-wavelength components of the Earth's static gravity field from GRACE (Gravity Recovery and Climate Experiment e.g., https://www. nasa.gov/mission_pages/Grace/index.html) and GOCE missions (Gravity field and steady-state Ocean Circulation Explorer, ESA, e.g., http://www.esa.int/ Our_Activities/Observing_the_Earth/GOCE). The "non-gravity" data to the RET 14 came from the Earth 2014 1 arcmin global topography model [Hirt and Rexer 2015] which incorporates the Bedmap 2 bedrock topography as a subset and the other data sets like topography from the Shuttle Radar Topography Mission (SRTM), Greenland Bedrock Topography and SRTM 30_PLUS bathymetry [for more details and references see Hirt et al. 2016].
EIGEN 6C4 [European Improved Gravity model of the Earth by New techniques, Foerste et al. 2014] is a combined global static gravitational solution, to degree and order 2190, including also a complete set of gradiometry data from the GOCE mission (in addition to the GRACE satellite-to-satellite tracking data and data from observations of other satellites).
Bedmap 2 [Fretwell et al. 2013] contains the bedrock elevation beneath the grounded ice sheet based on various types of ground, airborne and satellite data. The RES data are the most contributing. The bed topography in Bedmap 2 is given in a 1x1 km grid of heights of the bedrock above present day sea level (asl). The actual spatial resolution is, however, often lower than 1 km; the resolution and quality is still much worse for some areas without data or with sparse data and also the precision varies highly [see . Some data incorporated into the Bedmap 2 were taken from its predecessor known as Bedmap [Lythe et al. 2001, Figure 2b].
RET14 increases the resolution of EIGEN 6C4 and decreases the resolution of Bedmap 2; generally the spatial resolution of RET14 should be better than ~10 km over the whole Antarctica (excluding the polar gap and two bigger areas without data in Bedmap 2).
The success of our approach lies in the simultaneous use of all the gravity aspects and the bed topography and in the high resolution data now coming on the scene (for us today it is the RET 14). But even after that we can speak only about candidates for subglacial volcanoes, nothing more. The final decision about the existence of the volcanoes is expected from analyses by other geoscientists, namely after drilling the ice on the appropriate places; we provide geographic coordinates on the ice cover where to do it (Section 5.3).
The area of the Lake Vostok (LV) has been studied. Its bedrock topography, derived from Bedmap 2, is shown in Figure 1a. The LV is an elongated depression with two local minima. The hills in West and South directions (up and left in this Figure) continue to the lakes 90 East, Sovetskaya and others ("lake district"), via the Ridge B on the surface to the Western flanks of the GSM; to East and North, there is Wilkes Land (WL). The positions of our candidates for the subglacial volcanoes are encircled and numbered.

Reliability of the data
The data available were critically evaluated; we have a rich experience with various types of artefacts due to insufficient or missing data or their irregular, inhomogeneous distribution [see, e.g., Klokočník et al. 2016Klokočník et al. , 2017a], so we should be able to avoid misinterpretations also here in Antarctica.
The gravity data. The noise in RET 14 can be only roughly estimated to be better than 10 mGal nearly everywhere in Antarctica, based on information deduced from [Pavlis et al. 2012, Foerste et al. 2014, Fretwell et al. 2013, Hirt et al. 2016. These values are available only for the gravity anomalies [in the fundamental source Pavlis et al. [2012] they are called "commission error"]. But accounting for usually huge correlations among gravity anomalies inside the studied area and outside but not too far from it, we found about half values for the noise [Klokočník et al. 2017a]. These are our values N (as noise): N=5-10 mGal (more at the end of Section 4).
The bedrock topography data. How we can rely upon Figure 1a based on Bedmap 2? Can Examples 1 and 2 be real features, do they exist in nature or are they only artefacts generated by Bedmap 2 due to insufficient data? In Figure 1b we show the Bedmap 2 topography again (in black and white) and we add the tracks with RES, transformed from Figure 3 of Fretwell et al. (2013) into the geographic coordinates, and from Figure 2b of Lythe et al. [2001]. While LV is well covered by the data [with density of ground tracks from RES mostly in east-west oriented flights up to 7.5 km; see, e.g., Studinger et al. 2004], the tracks around the lake are not so dense (~ tens kilometres). In Figure 1b, the densely covered areas are shown in violet colour (the rectangles). They are there two large zones without data in Bedmap 2, called (by the authors of Bedmap 2) "poles of ignorance", PI; one is in Princess Elizabeth Land, see Fretwell et al. 2013, pp. 379 and 388. With our Example 2 we are close to it (Figures 1 a-c). We use RET 14 and it is influenced strongly by Bedmap 2. But they are some measurements exactly at or in vicinity of our Examples 1 and 2, namely RES data (see the tracks in Figure 1b). One older surface gravity traverse goes (by a good chance) directly over the subglacial mountains in the case of Example 2 in the older Bedmap (shown in Figure 1b as thick blue line); it is logical that such a data set was used also in Bedmap 2.
In the Bedmap 2 black and white topography, Figure 1b, we can see an evident trend to show details along the RES flight lines, changing often positive and negative heights quickly, so the lines with measurements look like "beads on a string" (Figure  1b). We can see a trend to smooth the bedrock topography among the tracks in places without data (Figure 1b, mostly its upper part, in west direction). But some topographic features are very pronounced so that they are not "averaged out" and that even sparser data can "catch" them.
A helpful fact is that Examples 1 and 2 are visible also by means of the gravity aspects from EIGEN 6C4 alone [Klokočník et al. 2016, Figure 17], but this test is strong only for Example 1 (more in Section 5). As we know, this gravity field model is completely independent of the both Bedmaps; that provides an independent and evidential check.
Finally a note to noise of Bedmap 2 heights. The estimated uncertainty in bed elevation grid is about 100 m at LV, but may reach 200 m in other parts we investigate [Fretwell et al. 2013].
We answer our question about reliability of the data: Example 1 and 2 might be real features, no artefacts.

Gravity symptoms of a single mountain from the gravity aspects
Sections 3 and 4 serve as a preparation for pres-entation of our results (Section 5). Here we focus on the particular object (a single mountain or a volcano) and list its typical gravity signal expressed in all the aspects. In each individual case of our testing of known volcanoes over the world, a single mountain, impact craters or other objects, the gravity aspects (listed above in Section 1) have a specific character and reach their local extrema at the objects and their vicinity. As a graphic representation to document our conclusions, we recommend to see Figures 3d-g and 4b-g of this paper and many figures in our previous papers [Klokočník et al. 2014[Klokočník et al. , 2016[Klokočník et al. , 2017b; they cannot be (due to space reasons) repeated here. A single mountain or volcano outside a mountain belt causes Δg, T zz , I 2 and I to have their local maxima; in the case of T zz the maximum at the top of the object is surrounded by a ring with negative values of T zz . For I 1 we observe a local minimum at its top. The quantities I 1 and I 2 have their extremes concentrated over the volcano's caldera. The volcano and its close vicinity exhibits a "missing θ". We compute θ for I < 0.3 and the localities with presence of flat causative bodies are indicated by one-side oriented values of θ; Figure 1b. Data gaps among the RES ground tracks (blue lines) in Bedmap 2 (and its predecessor Bedmap) plotted over the bedrock topography, taken from Bedmap 2 (in metres, asl) in geographic coordinates, near the Lake Vostok (LV). A gravity traverse from the older Bedmap [Lythe et al. 2001] is shown by thick blue line, see the upper right, just crossing our candidate for volcano Example 2. Our candidates for sub-glacial volcanoes are those two noticeable peaks aimed by red circles. Reproduced from [Fretwell et al. 2013] and adopted for our purpose. Figure 3 from [Fretwell et al. 2013] by anonymous reviewer (courtesy of ) to show positions of our candidates for subglacial volcanoes with respect to one of the "danger" (red) zones without data in Bedmap 2 [PI, Fretwell et al. 2013, pp. 379 and 388]. Scales of the x-and y-axes are in kilometres in a local system of Bedmap 2. for a 3D body like a mountain, the θ-values are missing.

The virtual deformations vd show a strong dilatation (red colour) at the caldera which is surrounded by a belt of compression (blue colour).
The gravity aspects and the bedrock topography must be accounted for together. If they all show what is expected for a volcano, then the tested mountain might be a volcano. But it is well known that the gravity information itself and alone does not provide unique results, so a special care is needed. The bedrock topography is needed, other data are welcome.
However, this way one cannot distinguish between a single "normal" mountain and a volcano. Thus, we modelled a volcano (Section 4); for its given shape, size and density, we computed the theoretical Δg M and T zz M for a reference point P above the top of the modelled object. Then we compared Δg M and T zz M with Δg O and T zz O of the actual volcano (the case of Mt Sidley, the largest volcano in Antarctica) and with the objects suspected to be a volcano, derived from the RET 14 (see below). The agreement of these values for Δg and T zz support our hypothesis about a volcano (Section 4).

Model of a volcano
This test should strengthen our predictions for the candidate volcanoes (objects) under the ice or under the sea level [for similar tests see Section 6 of Klokočník et al. 2016]. We developed our own simple model producing Δg and T zz values for a single mountain. We considered various input data about the single mountain (model of volcano), shape, size, density, etc. and computed the model values Δg M and T zz M . Then we computed Δg and T zz for the same object with RET 14, say Δg O and T zz O . Finally we compared both types of results.
The reference point P for our computations of Δg M and T zz M is put to the height h p vertically above the base of the object, always higher than all masses of the model (Figure 2). We consider the object as an isolated single peak (first important simplification); it has a circular base of diameter d and homogeneous density ρ (further simplifications). The volcano is modelled by a homogeneous truncated cone with a crater on its top. The ice cover is approximated by a layer with a uniform thickness h big enough that a truncated part would not affect the results. In the part of the volcano embedded in the ice, the density ρ is reduced by the density of ice to avoid a duplication of the ice contribution. A volcano body is divided into a set of mass elements, which are treated as mass points.
Their gravity attraction is computed numerically and summed over the whole body. We performed many tests with varying the input values h p , d, ρ, and h.
First we used an arbitrary known volcano in Antarctica as a test to compute its Δg M and T zz M and then we extrapolate to hypothetical volcano under the ice. We considered volcano Mt Sidley, the largest known non-active volcano in Antarctica (4.2 km above sea level) with the upper part protruding through the ice about 2 km high, in the place known as the Executive Committee Range (Marie Byrd Land); the central part of the caldera of Mt Sidley has these geodetic coordinates: φ = 77 0 03'41" S, λ =233 0 49' 52" E. We take this mountain as an isolated single peak. We estimate d = 60 km, height h v = 3.8 km over the surrounding area and density ρ= 3 g/cm 3 . The altitude h p = 4.8 km of the point P was used, i.e.  here see a note below.] A decrease of Δg M and T zz M with a distance from a centre of the volcano below the point P is sharper for T zz M than for Δg M ; this is expected for the second derivative; the half values of Δg M and T zz M in comparison with those for the top are found for d = 40 km and 20 km, respectively.
We can vary h p , d or ρ. Let us consider ρ i ≠ ρ. Then, under the ice (which has 1 g/cm 3 ), we can use the result for ρ = 3 g/cm 3 but with a multiple factor (ρ i -1)/(ρ -1), so for ρ i = 2.7g/cm 3 we get: (2.7 -1) / (3 -1) = 0.85, and Δg M = 75 mGal and T zz M = 80 E. This time the comparison with Δg O and T zz O from RET 14 comes out better. The previous test concerned Example 1; similarly successful was our test for Example 2.
These examples are not a prove but can be used as one of the arguments that the suspicious objects Examples 1 and 2 might be volcanoes. But they may be normal single mountains as well. We used the gravity symptoms (Section 2) and the bed topography (Section 3).
We care about a statistical significance of our results, too. The noise of RET 14 is N = 5-10 mGal [Section 2.2], following Sect. "Discrimination criterion" in Klokočník et al. [2017a]. The Signal (S) is taken from detailed listings belonging to the maps of the gravity disturbances over the objects 1 and 2. We seek for the S/N ratio; it should be > 3 to believe that we have statistically significant results. We find S/N = 7-14.

The results -Examples of possible subglacial volcanoes near the Lake Vostok
At the Lake Vostok, a seismic activity was recorded [Studinger et al. 2003b[Studinger et al. , 2004; thus we speculate that volcanoes might be expected not far away. We prefer a single mountain to a mountain in a mountain belt, because it can be determined with higher probability from existing gravito-topographic data. We could choose many other places in Antarctica and we certainly have a chance to find other candidates for subglacial volcanoes [Klokočník et al. 2016]. It is the highest non-active volcano known in Antarctica. The simple modelling of this volcano (see above) shows that our model can approximate well the real volcano. It means that the model works fairly and can be used in an extrapolation for the volcano candidates (next Section). The reader can follow the details below.

1. Test with known large volcano
All the symptoms (local maxima of Δg, T zz , I 2 , I; min. I 1 , no θ, positive vd at the volcano above its caldera, negative vd in a belt around the object) are fulfilled, as expected for each volcano or a "strong" single mountain. Here we show the bedrock topography from Bedmap 2 in Figures 3 a,  The topography (Figure 5a) reveals an isolated non-sharp-edged peak for Example 1, surrounded by a flatland and an oblong feature for Example 2 with a mountain belt nearby (Figure 5b).

Example 1.
There is a fair agreement between the actual Δg O and T zz O computed with RET 14 and the single mountain model for a volcano. The topography from Bedmap 2 confirms the existence of an isolated peak 1200 m above the relatively plain surrounding terrain (Figure 5a). Candidate 1 cannot be an artefact due to insufficient data in Bedmap 2 in this area, although the density of the relevant data here is not too high [see Figure 3 in Fretwell et al. 2013, our Figure 1b]; the reason confirming existence of Example 1 is that the same feature (a single massive mountain) computed with RET 14 (dependent of Bedmap 2) is shown (with a lower resolution) also by EIGEN 6C4 (independent of Bedmap 2), namely by T zz and vd [cf. Figure 17b in Klokočník et al. 2016]. Note also that the Bedmap 2 model used the data of the original Bedmap model [Lythe et al. 2001] and for this area there were some data already in this model [see Figure 2b in Lythe et       Conclusion. All pieces of circumstantial evidence about Example 1 lead to a conclusion that it is a candidate for a subglacial volcano (Figure 5a).        Figure 3 in Fretwell et al. 2013]; it is the LV area, but nearby in the north-west direction, an extensive "pole of ignorance" (PI, a large area still without any data in Bedmap 2) begins [Fretwell et al. 2013, pp. 379 and 388], see Figures 1a and c. One has to be very careful here to avoid a misinterpretation. The gravity data from EIGEN 6C4, which were helpful for Example 1, do not help here too much, because there is no single mountain surrounded by a plain area but a mountain belt; the resolution of EI-GEN 6C4 cannot distinguish the individual peaks and a smoothing takes place. Therefore, the older surface gravimetric data [a gravity traverse used in Bedmap, Lythe et al. 2001], just crossing our candidate for subglacial volcano Example 2, is very welcome (see Figure 1b right upper corner, thick blue line).

Geodetic coordinates and suggested names of the candidates for volcanoes
We provide the geographic coordinates, i.e. geodetic latitude φ and longitude λ (in WGS 84) for the top of these objects under the ice (Figures 5a,  b), according to the maps from Bedmap 2, together with the roughly estimated depths h of those tops beneath the ice layer given by the thickness of the local ice cover (derived from another source of Bedmap 2). These data might be used for prospective future drilling.

Conclusion
By using the bedrock topography Bedmap 2 and a high resolution static gravity field model for Antarctica RET 14, we computed the gravity anomalies/ disturbances, the Marussi tensor of the second derivatives of the disturbing potential, the gravity invariants and their specific ratio, the strike angle and the virtual deformations. We found by extensive testing that these gravity aspects show specific values for a massive single mountain. This knowledge was applied for our study area in Antarctica; not too far from the Lake Vostok we discovered at least two objects that might be subglacial volcanoes. We present all the arguments supporting this idea -the expected values of the gravity aspects using RET 14, topographic signal based on Bedmap 2 bedrock topography, and a good agreement between the gravity anomaly and the radial second derivative derived from RET 14 and from an independent simple volcano model. However, we are aware that the final decision about the respective volcanoes is on other specialists. We provide geodetic coordinates of the candidates for volcanoes (named Dana and Zuzana) and an estimate of the depth of their tops under the ice cover. Figure 5a. The bedrock topography from Bedmap 2 (m) of Example 1. One single peak with a relative height difference (with respect to surrounding landscape, mostly plain area) ~1000 m, the top has 1240 metres asl (sea level = black thick contour line with zero). A possible caldera below our resolution. Interval of contour lines 20 m, the lines in step 100 m are shown as thick curves; those below sea level are in blue. supported by the project LO 1506 (PUNTIS) from the Ministry of Education of the CR. We thank B. Bucha for the computations of the gravity aspects with EIGEN 6C4 and RET 14 and H. D. Pritchard for his help with the bedrock topography data.