On feasibility to detect volcanoes hidden under the ice of Antarctica via their “gravitational signal”

Many undiscovered volcanoes may be hidden under the thick layers of the ice of Antarctica. Hypothetic volcanoes were sought by means of the best present-day gravitational data (gravity field model EIGEN 6C4) and bedrock topography data (Bedmap2). A new previously unused method was tested. The analogy with the “gravitational signal” typical for volcanoes and other structures in other parts of the Earth is used. Various functions (not only ordinary gravity anomalies) of the disturbing geopotential were employed: Marussi tensor of the second derivatives, geopotential invariants, the strike angle and the virtual deformations. We attempted to discover if the best present-day gravitational and topographic data are of sufficient precision and resolution and how fast is the attenuation of the “gravitational signal” of a volcano with increasing depth under the ice. It is shown that there is no principal obstacle to detect volcanoes by our method. However, it appeared very quickly that our present-day attempts to discover such volcanoes could hardly be successful, mainly due to a low resolution of the existing gravity data and also due to a low resolution of the best bedrock topography of Antarctica currently available. Nevertheless, some examples of hypothetical volcanoes under the ice are given, but they are uncertain. However, the method, the main goal of this feasibility study, is ready and working.


Introduction
Volcanoes are nearly everywhere on the planet Earth (and on some other planets or their moons), so they are also in Antarctica (Figure 1).Many undiscovered volcanoes may be hidden under the thick layers of the ice (up to 4 km).They might be important for global warming, if some (huge) of them would become active.The discovery of two volcanoes active (smouldering) under the ice [e.g., Corr andVaughan 2008, Lough et al. 2013], new global and detailed gravitational field models with high resolution (like EGM 2008 [Pavlis et al. 2008a[Pavlis et al. , 2008b[Pavlis et al. , 2012] ] or EIGEN 6C4 [Förste et al. 2014]), derived from satellite as well as terrestrial data and the progress in mapping the topography of bedrock (Bedmap2 [Fretwell et al. 2013]) has been inspiring us to seek for hypothetic volcanoes hidden under the ice of Antarctica by using these data sources.An analogy with the "gravitational signal" known for volcanoes and other structures in other parts of the Earth was used.It is very important that (1) various functions (not only the ordinary gravity anomalies) of the disturbing geopotential (being represented by harmonic coefficients in the expansion of the potential to the spherical harmonic series, namely by EIGEN 6C4) were utilized (Section 2).It is also crucial that (2) we make use of the bedrock topography data (now Bedmap2, more in Section 3.2).These items (1) and ( 2) are two pillars of our method.
We discuss if the present-day gravitational and topographic data (Section 3) are of sufficient precision and resolution (Sections 3 and 5) and how fast is the attenuation of the "gravitational signal" of a volcano with its increasing depth under the ice (Section 6).It was easy to verify that the resolution is not sufficient now, but it is not a principal problem; it is anticipated that the new data will in the course of time provide updates of the Antarctic gravity data (Section 3).The attenuation with increasing depth of the ice is far from being critical (Section 6).
The methodology based on (1) and ( 2) is the main goal of this paper; the volcanoes are only the example.Our method is promising for a future successful search (with new forthcoming data of higher resolution) for subglacial volcanoes, having of course in hands also other than satellite data (Section 3).Our method is potentially valuable for all tectonic purposes and goes well beyond volcano spotting (partially already tested and proved successful for the areas elsewhere than in Antarctica, see Kalvoda et al. [2013] or Klokočník et al. [2014]).We need allies from geoscientists for this purpose for a further progress in our geoapplications.Detection of the subglacial lakes in Antarctica might be another goal (here we work only with the largerst one, the Lake Vostok, Section 7.3; a revised inventory of the lakes is in Siegert et al. [2005]).
Our present-day attempts (examples in this paper) to discover under-ice volcanoes can hardly be too successful, mainly due to a low resolution of the existing gravity data and also partly due to a low resolution of even the best bedrock topography of Antarctica now available (Bedmap2).But some tentative examples are given (Section 7.3).

Outline
The theory presented here comes from Pedersen and Rasmussen [1990], Beiki and Pedersen [2010] and from our recent papers Kalvoda et al. [2013] and Klokočník et al. [2014].A short outline of the theoryfor a convenience of the reader (non-specialist in physical geodesy) -is presented now.
The disturbing static gravitational potential outside the Earth masses in the spherical harmonic expansion is given by the well-known formula: where GM is a product of the universal gravitational constant and the mass of the Earth (also known from satellite analyses as the geocentric gravitational constant), r is the radial distance of an external point where T is computed, R is the radius of the Earth (which can be approximated by the semi-major axis of a reference ellipsoid), P l,m (sin {) are the Legendre associated functions, l and m are the degree and order of the harmonic expansion, ({, m) are the geocentric latitude and longitude, C' l,m and S l,m are the harmonic geopotential coefficients (Stokes parameters), fully normalized, C' l,m = C l,m − C el l,m , where C el l,m belongs to the reference ellipsoid.
The spherical approximation of the gravity anomaly Dg (on free air) is computed or one can use the gravity disturbance (it is the same as Equation 2 but without the second term, which is usually numerically small).The gravity anomalies/disturbances are computed from measurements by gravimeters or derived from measurements performed by means of satellite altimetry.
Gravity gradient tensor C (the Marussi tensor) is a tensor of the second derivatives of the disturbing potential T of the particular gravitational field model known to the maximum degree L max (see below some details about EGM 2008 and EIGEN 6C4).The second derivatives are expressed in the local north-oriented reference frame (x, y, z), where z has the geocentric radial direction, x points to the north, and y is directed to the west [Pedersen and Rasmussen 1990]: Outside of the source masses, C satisfies Laplace's differential equation, the trace of the Marussi tensor is zero, tensor C is symmetric, and it contains just five linearly independent components.These can conveniently be computed by means of the formulae in Hotine [1969].The tensor components are used in the local scales to identify and to map the geological contact information, be it the edges of the source targets or the structural/stratigraphic contact information.The horizontal components help to identify the shape and the geological setting of a target body.The quantity T zz is best suited for a target body detection; T zz helps to , , cos sin sin (3) define isopath/density relationships of a body mass with relation to its geological setting [e.g., Saad 2006].Under any coordinate transformation, C preserves three invariants and it is zero outside the masses of the studied body.Let us define (5) Pedersen and Rasmussen [1990] showed that the ratio I of I 1 and I 2 defined as (6) lies between zero and unity for any potential field.If the causative body is strictly 2D, then I is equal to zero.
The second order derivatives and the invariants provide evidence about details of near-surface (not deep) structures.The Marussi tensor was already used locally (areas of a few kilometres) for petroleum, metal, diamond, groundwater etc. explorations (e.g., Saad [2006], Mataragio and Kieley [2009], Murphy and Dickinson [2009] and many others).The full Marussi tensor is a rich source of information about the density anomalies providing useful details about objects located closely to the Earth's surface.This extra information can be used by tensor imaging techniques to enhance target anomalies, as tested for local features (economic minerals, oil and gas deposits, fault location, etc.), see, e.g., Murphy and Dickinson [2009], Saad [2006].
The strike angle i S (also known as strike lineaments or strike direction) is defined as (7) within a multiple of r/2.The strike angle indicates how gradiometer measurements rotate within the main directions of the underground structures.Provided that the ratio I (6) is small, the strike angle may indicate a dominant 2D structure.For more details see Beiki and Pedersen [2010] or Murphy and Dickinson [2009].
To define the new term "virtual deformation" [Kalvoda et al. 2013, Klokočník et al. 2014] (vd), an analogy with the tidal deformation was utilized; one can imagine directions of such a deformation due to "erosion" brought about solely by "gravity origin".If there would be a tidal potential T, then the horizontal shifts (deformations) would exist due to it and they could be expressed in north-south direction (latitude direction) as and in east-west direction (longitudinal direction) as where g is the gravity acceleration 9.81 m•s -2 , l S is the elastic coefficient (Shida number) expressing the elastic properties of the Earth as a planet (l S = 0.08), { and m are the geocentric coordinates (latitude and longitude) of a point P where we measure T; the potential T is expressed in [m 2 •s -2 ].In our case, T is represented by Equations ( 1), ( 8) and ( 9).This mechanism is applied to a standard Earth gravitational model, but the real values of the Shida parameters l S for the Earth's surface (for our purpose) are not known.
The apparatus of mechanics of continuum to derive the main directions of the tension is applied (see, e.g., Brdička et al. [2000]).The tensor of (small) deformation E is defined as a gradient of shift.It holds that The tensor E can be separated into two parts: (11) where e is the symmetrical tensor and X the anti-symmetrical tensor of deformation, respectively.The symmetrical tensor is: (9) a = ½ (D+ c) major semi-axis of ellipse of deformation b = ½ (D− c) minor semi-axis of ellipse of deformation a = ½ atan(c 2 /c 1 ) direction of main axis of deformation.
Note that different specializations have different terminology for the same or similar quantities.
To illustrate vd, the semi-axes of deformation ellipse a and b are computed together in their relative size.Values of l S are not known, and, therefore, only the main directions of the vd (and not its amplitudes) can be shown.The plotted quantities are a and b, expressed in our figures as crosses.The computation of all quantities defined above was organized by software based on Holmes et al. [2006], Sebera et al. [2013] and Bucha and Janák [2013].

Comments on the theory
A simplification -various properties of the aspects.To show how the various aspects (derivates, functions) of the disturbing geopotential, defined in the previous section, behave (decrease) with increasing distance from the source of causative body (a density anomaly), it is assumed that the source body is a mass point, with a selected value of GM, located on the Earth surface; all masses of the Earth are concentrated in the centre of the Earth.It is obvious from the definitions in Section 2.1 that in this simplified case, some aspects do not exist or equal zero.The remaining aspects are shown in Figure 2; here the simplification leads to the invariants I 1 = −3(GM) 2 r -6 and I 2 = −2 (GM) 3 r -9 .A diverse speed of decrease of the "gravitational signal" with increasing depth of the source mass body for various aspects can be seen.The slowest decrease is for the gravity anomaly (with r 2 ).The invariants decrease quickly (till with r 9 ).It tells us that aspects I 1 or I 2 best describe the density anomaly at the surface or in shallow depths under the surface while Dg and T zz are related to deeper struc-tures.While all these aspects have the same "common mother" (a recent gravity field model, Section 3.1.),they have different properties (behaviour, sensitivity).Thus, a set of the aspects from Section 2.1.provides much more information about the density anomalies than only the traditional Dg, or Dg together with T zz , can do.
Plotting figures -nonlinear scales.In our figures, strongly non-linear scales are used to emphasize various features which otherwise might remain hidden.The gravity anomalies and/or disturbances are in mil-liGals (mGal), the second order potential derivatives are in Eötvös; recall that 1 mGal = 10 -5 ms -2 , 1 E = 1 Eötvös = 10 -9 s -2 .The invariants I 1 and I 2 have units [s -4 ] and [s -6 ] and the ratio I is dimensionless.The strike angle i is expressed in degrees with respect to the local meridian and its demonstration in red means its direction to the west and in blue to the east of the meridian.These units are used in all presented figures.
Our figures for Antarctica are in stereographic polar projection with non-deformable parallel at 70°s outh (in this projection we obtained the topographic model Bedmap2).
Virtual deformations of the ellipse of deformation (vd) are geometrically expressed by its dilatation (shown in red colour) or compression (blue).The virtual dilatation indicates uplifted regions at the geoid, whose mass has (owing or according to the pattern of values of the gravitational potential) a tendency to disintegrate.On the contrary, the virtual compression indicates lowered zones at the geoid.Natural processes, which are the cause of these features near the surface of the geoid, are certainly very diverse as a consequence of regionally heterogeneous integration of morphotectonic, erosion-denudation processes, and others.
A systematic screening of the whole planet Earth was performed.The results of screening were represented by means of selected examples from regions of various planation surfaces, high mountain ranges, collision zones of oceanic and continental lithospheric plates, regional fault zones, volcanic chains and large impact craters, in Kalvoda et al. [2013], Klokočník et al. [2010Klokočník et al. [ , 2014]], or Klokočník and Kostelecký [2014] (www.asu.cas.cz/~jklokocn).Here we present typical examples of the "gravitational signal", primarily for the volcanoes (Section 4).

Gravitational field models with high resolution EGM 2008 and EIGEN 6C4.
The input data to our analysis are the harmonic geopotential coefficients (also known as Stokes parameters) in the spherical harmonic expansion of the disturbing gravitational potential T, Figure 2.An example of decrease in the values of the aspects of the disturbing geopotential with increasing distance (depth) from their source.On the x axis there are kilometres, on the y axis there is an arbitrary quantity (for a simple intercomparison).This is an illustrative case for a mass point, representing a density anomaly, with a random value of GM.
namely numbers C lm , S lm in Equation ( 1), Section 2.1.A set of these coefficients defines a gravitational field model.Today, global gravitational field models of the Earth, based on variety of satellite and terrestrial data, have a high resolution.The Earth Gravitational Model 2008 (EGM 2008[Pavlis et al. 2008a, 2008b, 2012]) or European Improved Gravity model of the Earth by New techniques (EIGEN 6C4 [Förste et al. 2014]) are both expanded to degree and order 2190 in spherical harmonics.This corresponds to a resolution of 5×5 arc minutes, which is ~9 km half-wavelength on the Earth's surface (details below).In this paper, we used mostly EIGEN 64C -always for Antarctica.For more details on precision, see Pavlis et al. [2012].With EIGEN 6C4, the "gravitational signal" is more precise for the majority of parts in Antarctica than with EGM 2008 due to the GOCE data in EIGEN 6C4.
Unfortunately, Antarctica is an exception from that limit of 9 km; the models have here a lower resolution, because no terrestrial gravity data are included (although some are available in few areas of Antarctica).It means that for Antarctica, only satellite gravity data are available.They are from GRACE [e.g., Tapley et al. 2004] and GOCE measurements [e.g., Floberghagen et al. 2011].GRACE is an acronym for Gravity Recovery And Climate Experiment, GOCE for Gravity field and steady-state Ocean Circulation Explorer, ESA's gravity mission, for the first time with a gradiometer on board (measuring directly the components of the Marussi tensor, Equation ( 3)).For more details about the accuracy and resolution of the "gravitational signal" over Antarctica see below and in Section 5.
Polar gap.The satellite data in the gravitational models are not available over the whole globe but have small polar gaps due to the fact that satellite orbits are inclined to the Earth's equator with an angle different from 90°.GRACE has an orbit inclination of 89° and GOCE of 96.7°, so that the polar gap without data is as large as 1° and ~7° around the poles for the GRACE and GOCE, respectively.Because the GRACE data have a lower resolution than the GOCE data, then, going from the ocean to Antarctica, first we lose the precision and resolution having only GOCE+GRACE data (and no terrestrial data) and then, very close to the poles, again we lose precision because we have only the GRACE data (see also Figure 9a).The narrow polar gap 1° around the southern pole is without any data in EIGEN 6C4 and the aspects (2)-( 7) and vd should be ignored here.
Theoretical resolution.The expansion to degree and order 2190 means a resolution about 10 km on the ground (see above).More precisely: a representation of a function in spherical harmonics has always some practical limit in degrees L max < ∞ and it corresponds to a spatial resolution at the Earth surface.A usual simple estimation of the smallest representable feature of the gravity field or the shortest half-wavelength L half (as a distance on the sphere) that can be resolved with all C lm , S lm to L max , is or equivalently, taking the circumference 2rR=40,050 km: For L max =300 and 2190, we have L half =67 and 9 km, respectively, which will be used later (Section 5).The former limit corresponds to the resolution of GOCE data separately and the latter characterizes the combined gravity field models EGM 2008 or EIGEN 6C4.
Good to note that the theoretical, "formal" or "geometrical" criterion ( 13) is just first "approximation to the truth".In reality, the resolution depends on the quality of data, their distribution over the globe, and on data processing.It should not be so low (67 km only) for Antarctica; more to this topic see in Section 5, Figures 9a-c.
The resolution of the gravity data today is obviously the main practical problem for our search for volcanoes hidden under the ice of Antarctica.One cannot find any individual volcano (even such huge as Mauna Loa/Kea) in this way, but one can combine the "gravitational signals" (Section 2.1) with the bedrock topography (Section 3.2).Recall that we have not only the gravity anomalies available for our search, but we work with other aspects (functions) of the disturbing geopotential T, as defined in Section 3.1.They are helpful because of their different sensitivity of resolving power (higher than that of the gravity anomalies) to the density contrasts located shallowly beneath the surface (remember Figure 1).But they still do not have the required resolution, it is obvious.Tests of resolution follow in Section 5.

Bedrock topography data
The topography of the ground under the ice (bedrock, base of the ice sheets) in Antarctica, known as Bedmap2, was compiled from 25 million of measurements of various kinds with a resolution reaching 1×1 km [Fretwell et al. 2013], but not everywhere.Note the ice thickness in Bedmap2 used to make the bedrock topography across most of the ice covered regions was gridded at 5 km; this grid was then re-sampled at 1 km resolution.As the original gridding was done on a 5 km grid, this should be considered as the resolution of Bedmap2 in regions of extensive ice cover.
The data ("bedmap2_bed") contains the bedrock elevation beneath the grounded ice sheet of Antarctica in the grid 1×1 km in WGS 84 in the polar stereographic projection; the data were downloaded from [Fretwell et al. 2013].Errors in ice sheet bed elevation are ~60 m in regions where there is data, rising to 150 to 200 m in regions further from observations, and may be >1000 m in regions >50 km from any direct measurements.This is a bit risky for user because he/she does not know exactly the quality of data in the particular locality.
The subglacial surface derived from the Bedmap2 data is used here together with the gravitational data mentioned above in Section 3.1 (but will not be com-bined with them to create for example Bouguer or isostatic gravity anomalies, in a contrast to, e.g., Bonvalot et al. [2012]).

Other data
A new degree-2190 gravity field model Sat-GravRET2014 (shortly RET14) in terms of C lm , S lm , valid or meaningful from physical viewpoint only for Antarctica, has been published very recently [Hirt et al. 2016]; it combines EIGEN 6C4 and Bedmap2.It is obvious that RET14 increases the resolution of EIGEN 6C4 and decreases the resolution of Bedmap2; generally its resolution should be ~10 km (like for the EIGEN 6C4 it-  self outside Antarctica).We computed all the quantities defined in the section "Theory" for the same regions of Antarctica as with the EIGEN 6C4 but we used these results only to a limited extent.We guess that for our purpose it is better to keep the gravity (EIGEN 6C4) and the topography (Bedmap2) separate and mutually independent.
Other important data sets have appeared recently, but not in terms of C lm , S lm , thus not directly applicable with our theory.Operation IceBridge was a 2009-2016 NASA mission aiming to monitor changes in polar ice from a fixed-wing aircraft [e.g., Studinger et al. 2010].It is a temporary replacement for the ICESat until ICE-Sat-2 launch in 2017.The new Antarctic gravity anomaly grid [Scheinert et al. 2016] incorporates all available gravity data over Antarctica over the last three decades.The PolarGap project [Forsberg et al. 2016] nearly filled the GOCE polar gap (Section 3.1) with the airborne gravity, magnetic and radar data, reaching precision 1-2 mGal rms in Dg.

Typical gravitational signal of various geological structures, namely volcanoes
The quantities ( 2)-( 7) and vd were computed and few examples were plotted to learn how a typical signal from volcanoes (and other structures) looks like.This knowledge is later (below) extrapolated for Antarctica.The questions are: what are the characteristics of the gravity field around a volcano?How to distinguish hypothetical volcano from an "ordinary" mountain?
Figures 3a-c show Mount Fuji area in Japan, the radial derivate T zz (Equation 3) and the invariants I 1 , I 2 (4, 5).Fuji has a large positive value of T zz .The invariants are still more sensitive than T zz and exhibit strong variations of the signal around Fuji and surrounding mountains.Figure 3d   in the Pacific ocean near Japan; note the strong positive vd inside the volcanoes surrounded by a "moat" of negative values, all extremely locally concentrated.Here, under water, the signal is certainly partly attenuated (see Section 6 for more details concerning the signal under the ice).
Figures 4a-e show the Big Island of Hawaii; one can see details for the individual volcanoes and lava flows using T zz and vd, while the gravity anomalies reveal only a basic shape; the dilatation in vd over the island prevails; a belt of the depression of vd is formed around the island.For the oceanic examples (Hawaii and in the Pacific) the signatures most likely reflect isostatic compensation of the volcanic loads.
Figures 5a-b show Popocatepetl and Iztaccihuatl near Mexico City in terms of T zz and I 2 .It is typical that large volcanoes have the strongest gravity signal (positive values) of all mountains in the close vicinity.But well visible are also "valleys" or "moats" around the volcanoes (negative values of T zz ), which in some cases can be (together with stripe-like structures nearby) an artifact of gridding; if the causative body is small, approaching in size the resolution limit of the used expansion (here L max = 2190 or 9 km), then a "bull-eye" effect and the "stripes" may appear.The examples are in Figures 5c-d.In 5c, we have Mt.Olymp on the planet Mars, the largest (shield) volcano in the solar system.But the current best available gravitational field model available for Mars [Marty et al. 2009] is expanded "only" to degree and order 95 (yielding a resolution of ~200 km), so even such a large object is "small" in comparison with the resolution limit and the artifacts immediately appear (the blue zones do not exist in reality).Figure 5d shows two smaller volcanoes on the Egyptian-Libyan border, with the same artifact effect (using the resolution of the full EIGEN 6C4, which means ~9 km on the Earth surface).
Figures 6a-c are devoted to the gravitational signal in Israel and surrounding countries, confirming, according to geophysicist L. Eppelbaum (Tel Aviv University), the well-known structures in this area (more in Klokočník et al. [2014]).It provides one of the testbeds for extrapolations to other places.
The large impact crater Popigai, Siberia (possibly with its smaller companions in SE-NW direction) is in Figures 7a-c  aly.(We note that Popigai is a place of large diamond deposits, created by the impact).
It is demonstrated in this section that a large volcano typically generates strong gravity signal as high positive values of Dg and T z , together with extreme values of the invariants I 1 , I 2 and positive vd "in the volcano" (dilatation), surrounded by a belt-like feature of negative vd (compression).But mountain belts without volcanoes can also generate a large gravity signal (not shown here).It is not easy to distinguish volcanoes, it is not unambiguous, but the volcanoes have very strong, concentrated signal and their I 1 , I 2 and vd show typical sharp extremes.The extrapolation of this knowledge to Antarctica it is not so straightforward (more in Sections 5 and 6).If it were easy, it certainly would be already solved and published.In Antarctica, we have to count with a lower resolution of EIGEN 6C4.

Problem with the resolution of the gravitational data
The problem of the insufficient resolution of EIGEN 6C4 (and similar global gravity field models), with L half ~70 km for Antarctica, Section 3.1, looks like a serious obstacle precluding a possibility to discover any, even a very large volcano under the ice of Antarctica.Let us investigate the resolution of the recent EIGEN 6C4 (or EGM 2008) more closely.
We performed tests to model the resolution corresponding only to satellite data over Antarctica (max.degree 300) instead of the full resolution (2190).EIGEN 6C4 was cut at degree and order 300, then the functions (2)-( 7) were computed and plotted again for a selected volcano or their group; one example is shown for the Popocatepetl area, namely T zz (Figures 8a, to be compared to Figure 5a) and two examples (Dg and T zz ) for Hawaii (Figures 8b-c to be compared to Figures 4a-d peared (Figure 8a).For the large Hawaiian shield volcanoes one can see only their main feature (a positive anomaly of circular shape).This "smoothing effect" will negatively affect our activities described in Section 7.3.Note that for maximum degree 300, there is again the aliasing gridding effect (as in Section 4), clearly seen in Figure 8b in a form of stripes arranged roughly along the islands of Hawaii.
It is evident (Equation 13) that with the resolution up to degree ~300 only, one cannot discover any volcano under the ice of Antarctica.But there is an additional information to the simple criterion in Equation ( 13) from Section 3.1; this may shift the limit from (13) to a "better resolution".First, one can combine the gravitational data with the topography from Bedmap2 with higher resolution [e.g.Bonvalot et al. 2012].But this is not our way, because we wish to keep the topography as an independent source.
Second, a better insight into the quality of satellite measurements over the polar regions (excluding the polar gaps themselves which are free of data) yield Figures 9a-b.In Figure 9a, the reader can see a distribution of the ground tracks of GRACE and GOCE satellite orbits over Antarctica and in the large range of southern latitudes outside Antarctica.A significant densification of the tracks in the direction from smaller latitudes towards the pole can be seen.Therefore, the gravity signal based on such data should be better sampled going to the pole, compared with that over the equatorial areas, reaching the maximum of the coverage density  at the turning points of the satellite orbit (corresponding to the orbit inclination of the respective satellite) and then reaching immediately zero (the polar gap for the given orbit).
Figure 9a shows the densification of the satellite ground tracks near the turning points of GRACE (red) and GOCE (blue).Under typical orbital conditions, the quantitative increase of the number of the ground tracks per unit area is shown in Figure 9b.In this figure, the surface density was normalized by the density at the equator (where its value is set equal to one).The inclination of GRACE orbit is ~89°.The GRACE observations (in red) cover the Earth surface almost over the South Pole, the coverage leaves a 1° gap around the pole, Figure 9a.The density of the ground coverage is shown to increase towards the pole; Figure 9b demonstrates that in the latitude band above the turning point there are approximately 100 times more GRACE observations per unit area than over the equator.At latitudes 70°-80°, there are ~6 times more observations,  Here the surface density is defined as the number of observations per unit area and it is displayed as function of latitude (here the density is normalized by its value on the equator).The y-axis is shown linear (upper plot) or logarithmic (lower plot).(c) Geoid height error (in meters), as a function of geographic location, of a satelliteonly gravity field solution ASU-CHAMP-0309 computed from 7 years of CHAMP GPS data (from Bezděk et al. [2014]).
thus the observations should produce √6 ≈ 2.4 times better results in precision of the estimated spherical harmonic coefficients.In order to extract the maximum information from observations, usually the gravity solutions are computed to higher spherical harmonic degrees.Then, at some harmonic degree L max the observation noise overweights the signal and the spatial resolution of the model may be characterized through Equation ( 13).If the number of observations is globally increased, then a better signal-to-noise ratio implicates higher limit L max , thus giving better spatial resolution, too.This idea can be extended to regions with more precise observations (or more dense observations), where the fitted model produces more precise coefficients for higher harmonics, and thus here the spatial resolution may be (somewhat) higher than predicted by the rule ( 13) giving an overall characteristic of the model.The inclination of GOCE orbit was ~96.7°, the consequence of which is a larger polar gap of 6.7° (see Figure 9a, curves in blue).The turning point of GOCE is located closer to −80°, where there are ~9 times more observations per unit area than at the equator, thus giving an improvement factor ~3.
The densification of the ground track coverage near the orbital turning points is further illustrated on an example derived from real satellite data.In Figure 9c, the geographic distribution of the error in the gravity field model computed from the GPS positions of CHAMP is shown (for details, see Bezděk et al. [2014]).It is clear that precision and resolution of the model is increasing with latitude towards the poles.For models based on satellite observations, this latitude-band pattern is typical [cf., Tapley et al. 2005, Pavlis et al. 2012].
In near future, a substantial improvement of the bedrock topography data and maps, further extension of seismic data, ground (here ice) penetrating data as well as local (airborne) and satellite gravity data is expected.Earth's gravity field tailored missions like GRACE Follow-On (e.g., http://gracefo.jpl.-nasa.gov/mission/; https://www.ae.utexas.edu/news/features/grace-follow-on) and GOCE FO (http://meetingorganizer.copernicus.org/EGU2015/orals/17191,etc.) are in preparation; all this is leading to a higher precision and resolution in topography and gravity (also) over Antarctica.

Modelling the gravitational signal of a volcano under the ice
In Antarctica the gravity signature of sub-ice topography will be attenuated due to the decreased density contrast between the rock and surrounding medium.
To estimate the damping due to the ice layer of vari-ous thickness, the volcanoes are modelled by a homogeneous truncated cone, with a "crater" of diameter 0.5 km on its top.
The ice cover is approximated by a layer with a uniform thickness h, big enough that the truncated parts would not affect the results.The ice height h will vary between 0 and 2 km in our computations.
The slopes of the volcano model are taken as 40 and 9 degrees for a stratovolcano and a shield volcano, respectively.The base of the cone depends on the volcano height and the slope.A constant density of 2.6 tons per m 3 is used, but in the part of the volcano embedded in ice, it is reduced by the density of ice to avoid duplication of the ice contribution.The effect of the model is obtained by numerical integration.
The anomalies are calculated for a set of points distributed at an altitude of h P = 2.0 km over the volcano basis, which corresponds to maximum of either volcano or ice height (Table 1).Other altitudes and other density contrasts (the volcano vs. the ice) have also been tested (not in Table 1).
The results summarized in the table show peak values of Dg and T zz , while both quantities decrease quickly with the distance from the volcano.As the ice layer is homogeneous, it only yields a constant shift in the Dg over a large area around the volcano ("ice cap" in the table) and a change of T zz is quite negligible.Thus the only effect of the ice cover arises from the smaller relative density of the submersed part of the volcano.
The modelling has shown that there is no essential issue with detecting volcanoes covered partly or even fully by ice due to the ice damping of Dg and T zz .The attenuation does not exceed 40% of the undisturbed value (when h = 0).This type of attenuation is approximately equivalent to that seen in the ocean, where submerged volcanic seamounts are easily observed.

General view over the whole Antarctica
The quantities ( 2)-( 7) and vd were computed and plotted for the whole planet, in 5×5' window, including Antarctica.Examples in a form of "general views" over the whole Antarctica are presented in this section, Figure 11a for Dg, Figure 11b for T zz , Figure 11c for I 1 , Figure 11d for I, Figure 11e for i, and Figure 11f for vd.The topography of bedrock from Bedmap2 is added as Figure 11g (including the names of the main topographic structures).Note (Figure 11a) the gravity anomaly to be low at the Ross sea, at Victoria and George V Land; next, the anomaly is crossing the continent over the South Pole from SE to NW, another is at Wilkes Land including a semi-circular structure, etc.There are two large and long positive gravity anomaly belts on the north side of Antarctica (separated by a long belt of negative anomaly), on the west at the Antarctic Peninsula and near the Gamburtsev Moun-tains.These features are very well supported and more details are provided with sharply defined boundaries by the values of T zz (Figure 11b).
ANTARCTICA UNDER ICE BY GRAVITY SIGNAL ice height (km)   Remarkable is a step decrease of resolution when going from the oceans around Antarctica (on oceans mostly the data from satellite altimetry were used in EIGEN 6C4) to the land of Antarctica (see explanation above).The second decrease in the resolution at ~83°b etween the zones with the GOCE+GRACE and GRACE-only data in EIGEN 6C4 (explained above) is not much noticeable.
The ratio I is in Figure 11d; the strike angle i (Figure 11e) for I<0.3 indicates, according to the theory (Section 2.1), the main directions of the underground structures and possibility of existence of a 2-D (read: flat) structures.It might be interesting for another study (by specialists, not by us) to correlate the blue or red zones of i to the places with large glaciers; for example two large blue zones correspond to Lambert Glacier, a major glacier in East Antarctica (symbol LG), and to Dehmann Glacier (DG), both shown by arrows in Figure 11e.
The values of vd over Antarctica are in Figure 11f using a reduced resolution to be able to plot this comprehensive figure (for more details, see below the figures with vd depicted for individual segments and localities).The areas with dilatation are in red colour, those of compression in blue.The strongest stresses appear to follow the continental-ocean transition, and elevated topography of the Transantarctic Mountains.They are  located also in a belt going along longitude 160° roughly from the south to the vicinity of the pole, at the Antarctic Peninsula (west) and near the South Pole (where, however, the data have the lowest resolution).
A general view on bedrock topography over the whole Antarctica (Figure 11g) is a fascinating result but here it does not bring anything new to us.Much more details for the selected areas are needed.Thus, the full resolution of the Bedmap2 maps was used; it should be 5×5 km (Section 3.1), see the following sections for various details.

Selected segments and places of Antarctica with known volcanoes
Figures 12a-d (for more examples see www.asu.cas.cz/~jklokocn) are not our final goal, of course.According to the distribution of the main known volca-noes of Antarctica (Figure 1) and accounting for our interests to try a search for volcanoes still hidden under the ice, few segments and places of Antarctica for a more detailed inspection were selected; plots of the functions of the disturbing potential represented by EIGEN 6C4 and by the Bedmap2-derived topography of the bedrock are presented below.
Although Antarctic volcanoes are technically stratovolcanoes, the frequent development of lava-fed deltas has resulted in the volcano profiles which are "flat", normally associated with lava-dominated shied volcanoes.
Segment 1: It includes the Ross shelf glacier and Ross Island (with volcanoes Mt.Erebus, active volcano, Mt.Bird and Terra Nova), Victoria Land and a part of Transantarctic Mountains (with volcanoes like Roy.Soc., Mt.Morning or Mt.Discovery).Figures 12a-d ANTARCTICA UNDER ICE BY GRAVITY SIGNAL  show T zz , I 2 and the topography from Bedmap2 for this segment.
In Section 4, it is demonstrated that volcanoes are located usually at the maxima of T zz and at the extreme values of the invariants.But this is not an unambiguous guide to find a volcano under the ice because in the areas with lower resolution of the gravitational signal, e.g., just in Antarctica, there is an additional complication.It was called a "smoothing effect" in Section 5 and the comparison of Figure 5a with Figure 8a or Figure 4d with Figure 8c implies how it works.The details of the gravitational signal disappear completely and/or the surfaces showing the gravitational signal are smoothed, shifting a location of a possible volcano from the maximum of T zz to the slope between the maximum and minimum of it, from a point value to a broader region.Moreover, Dg or T zz due to a nearby mountain belt or another geological feature can dominate the signal of the volcano itself.This happened for example for the Ross Island with T zz (Figure 12a) thanks to a close and large pos-itive T zz of a long mountain belt.The RAT14 (Section 3.3) provides here satisfactory results (Figure 12c).
From all known volcanoes in Segment 1 (according to Figure 1), 54% or 64% are located inside zones with positive values of Dg or T zz , respectively, 82% or 64% where I 1 or I 2 have extreme values, 64% in the case of vd, and 100% for the local peaks shown by the topography with Bedmap2 (visible volcanoes can be confronted with © Google Earth satellite images).
Segment 2: Mt.Sidley and the Executive Committee Range (ECR).There is a volcanic mountain range with the largest volcano in Antarctica (Mt.Sidley, mostly covered by snow, a well visible caldera with a diameter of ~5 km), together with Mt.Waesche, Mt.Hampton and others, creating a long range (see Figure 13a).The ECR is about 100 km long, in the NS direction.Figures 13b-g show examples of quantities ( 2)-( 7), vd and the topography from Bedmap2.The gravitational signal does not exhibit anything like a range (belt), but all volcanoes of the ECR are located on places where Dg and T zz are positive.The resolution of Bedmap2 on Mt.Sidley was tested, see the details in Figure 13f.It looks like that the promised resolution up to 5 km or even 1 km is not achieved here (Section 3.2) because we do not see the caldera, which has ~5 km in diameter (Figure 13g).But this is not a critique of Bedmap2, which we found to be a very useful tool; for our purpose, however, its resolution is not sufficient.

Search for volcanoes hidden under the ice of Antarctica
Segment 3: Gamburtsev Mountains.Here there are large mountains (like the Alps) with peaks, river valleys, etc., all "quickly frozen" and preserved "in a good shape" probably last ~14 million years in a "deep freezer", under an up-to-3.5 km thick ice layer [e.g., Bo et al. 2009].Figures 14a-f show the topography of this area from Bedmap2, then T zz , I 1 , i, I, and vd.We have encircled extreme values of that topography and places with no signal in i (for I <0.3), meaning that there should not be a flat body beneath the surface.Also a good sign.
The values of T zz , T yy and I delineate locations of extreme values from the north, east, west and south; only the ratio I is shown here, see Figure 14e; its large value (I > 0.8) in the encircled area cannot be hiding a flat body.The dilatation in vd in this place is also evident (Figure 14f ).This is a place with higher probability of an appearance of a volcano but preliminary geological interpretation [Ferraccioli et al. 2011] does not indicate any volcano.
Segment 4. This is a part of Queen (Dronning) Maud Land.The Google Earth photos show nunataks and mountain belts (rocky peaks, some exceeding 3600 m above sea level, piercing the ice cap) and Bedmap2 discloses a very rich bedrock topography (Figure 15a).We selected the area between roughly { = 70 and 75°S, m = 0 and 20°E.The area has only recently been explored as for the gravity data [e.g., Riedel et al. 2012].
Figures 15b-e show the relevant T zz , I 2 , i and vd with coincident extremes.Then details from Bedmap2 for selected features (see numbers 1-6 in Figure 16a) ANTARCTICA UNDER ICE BY GRAVITY SIGNAL follow (Figures 16b-d).
Here three examples of the detailed topography are shown and they suggest, in a correlation with the gravitational signal in Figures 15b-e, where possible volcanoes might be located.The "smoothing effect" (Section 5) may, however, push the maxima of the "gravitational signal" in such a way that volcanoes should be expected rather or also at the slopes of the signal between their maximum and minimum values than in the maxima or extremes directly.
Segment 5: Around Lake Vostok.Topography is shown in Figure 17a where one can see the lake as a large depression.It corresponds with the gravity signal; here T zz is shown (Figure 17b) and it has there a negative value, as expected.According to T zz the lake might continue in the south-east direction.Important is the isolated single "round" mountain in the central part of Figure 17a.This feature, resembling the shape of a volcano, is supported by a positive T zz , which is also expected for a volcano (17b, encircled on both figures) and is known from other parts of the Earth.May be another candidate for a volcano under the ice was discovered.But we do not know how many data in the Bedmap2 were available here, how they are precise, and how they were distributed geographically in this locality; we hope that the conic feature in Figure 17a is not just an artifact of Bedmap2.
There are many possibilities where volcanoes under the ice of Antarctica can be located, and as readers know from the previous text, it is mainly an insufficient resolution of the gravitational data what are now available what precludes their reliable identification.

Conclusions
This is a feasibility study, nothing more.It concerns methodology and examples rather than the particular results now attainable with existing data.
Theory, data and computations were presented in Sections 2, 3, and 7, respectively, in a form of figures for various aspects (functions) of the disturbing gravitational potential, represented by the Earth gravitational field model EIGEN 6C4 (probably the best combined model now (2016) available, with a full set of GOCE gradiometric data; Förste et al. [2014]).It is very important that (1) the various functions (not only the ordinary gravity anomalies) of the disturbing geopotential are used, and that (2) we make use of the bedrock topography data (Bedmap2) together with the gravitational data.The items (1) and ( 2) are two pillars of our technique.
Attenuation (damping) of the "gravitational signal" of a test (model) volcano under the ice layers is not critical (Section 6); the attenuation with increasing depth of the model under ice is below 40% of the undisturbed anomaly for both the gravity anomaly and the second radial derivative (the attenuation is a function of the density contrast (volcano vs ice), decreasing with increasing contrast).It means that enough of the signal generated by a large volcano under the ice remains to be detectable.Our tests, therefore, show that the attenuation does not pose any principal obstacle for a successful search of volcanoes under the ice of Antarctica.
Resolution of the recent gravity field models (like EGM 2008 or EIGEN 6C4) is about 10 km on the ground, but for Antarctica, where only the satellite data  from GRACE and GOCE are available, the resolution is lower (between 10 and 67 km, Section 5).The resolution of the under-ice topography with Bedmap2 is much better (it should reach 1 km, but not everywhere, 5 km is expected realistically) than that of the gravity data, but also not sufficient yet for our purpose.This does not, however, mean any principal obstacle to solve the given task in the future, but better data are needed.They may come from satellites GRACE and GOCE Follow-On missions and from new spaceborne laser altimeters as well as from terrestrial and airborne data (seismic measurements, ice penetrating radar, airborne gravimeters…).For example Forsberg et al. [2016] showed filling the polar gaps by airborne data.It looks like that next few years can bring considerable progress in our knowledge of the gravity field and bedrock topography of Antarctica.
A temptation to predict locations of possible under-ice volcanoes, using the gravitational and topographic data now available, was very strong, although a great success with the input data now available cannot   be expected.In Section 7, the "gravitational signal" and the bedrock topography for selected segments of Antarctica were presented, two examples for localities with known volcanoes and four predictions (extrapolations) of hypothetical volcanoes still masked under the ice (Figures 14-17).
Our method is potentially valuable for all tectonic purposes and goes well beyond volcano spotting in Antarctica.(As for the topography outside Antarctica, they are available such sources like ETOPO 1 or ASTER GDEM.)Detection of the subglacial lakes might be another application of our approach.
S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S (12) and the parameters of deformation are: D = e 11 + e 22 total dilatation c 1 = e 11 − e 22 pure cut c 2 = 2e 12 technical cut c

Figure 3 .
Figure 3.The gravitational signal of a volcano -type example.(a) The radial derivative T zz -note the large scale (more than 300 E); (b,c) the invariants I 1 , I 2 around Fuji; (d) the virtual deformations over Japan, at a close trench in the Pacific ocean and undersea volcanoes; computed with EGM 2008 till degree and order 2190.Reproduced from Klokočník et al. [2014].
Figure 4. Hawaii: (a) T zz for all the islands; (b) the Big Island, taken from © Google Earth; then for the Big Island alone (c) Dg, (d) T zz and (e) vd.
Figure 5. (a-b) The largest stratovolcanoes Popocatepetl (P) and Iztaccihuatl (I) near Mexico City; (a) T zz and (b) I 2 .(c) The values of T zz for Mt.Olymp on Mars.An excellent example of the gridding aliasing effect (bull eye and stripes), partly visible on Figure 5a as well.(d) T zz for two volcanoes (Uweinat and Arkenu) on the Egyptian-Libyan border; another example of the artifacts (bull eye).
Figure 6.(a) Dg, (b) I 1 and (c) vd for Israel confirming the structures well-known to geologists and geophysicists.Reproduced from Klokočník et al. [2014].

Figure 7 .
Figure 7. (a) Dg, (b) T z and (c) vd the for impact crater Popigai, Siberia (diameter of the main crater on the NW side is about 100 km).Reproduced partly (panels (a) and (b)) from Klokočník et al. [2010].

Figures 8
Figures 8 (same column, top to bottom).The truncation of EIGEN 6C4 at degree and order 300 to simulate a model with satellite (GOCE) data only.The values of T zz around Popocatepetl (a) and at Hawaiian Islands (b) and Dg for the Big Island of Hawaii (c).To be compared to Figures 5a and 4a,d.

Figure 9
Figure 9 (same, column, top to bottom).(a) Ground tracks of satellite missions GRACE (red) and GOCE (blue) near the South Pole.The latitude turning points result from orbit inclination, which is 89° for GRACE and 96.7° for GOCE.The ground track points of both missions are displayed for 5 days in December 2010.(b) Surface density of observation points for GRACE (red) and GOCE (blue).Here the surface density is defined as the number of observations per unit area and it is displayed as function of latitude (here the density is normalized by its value on the equator).The y-axis is shown linear (upper plot) or logarithmic (lower plot).(c) Geoid height error (in meters), as a function of geographic location, of a satelliteonly gravity field solution ASU-CHAMP-0309 computed from 7 years of CHAMP GPS data (fromBezděk et al. [2014]).

Figure 10 .
Figure 10.A model of volcano for investigation of the attenuation of the gravitational signal (in Dg and T zz ) under a layer of ice, as observed from the point P.

Figure 11 .
Figure 11.(a-b) The values of Dg and T zz over Antarctica with complete EIGEN 6C4.White dots are for the largest known volcanoes in Antarctica, mostly near the sea shore or at islands near the continent.An extra map of known volcanoes was in Figure 1.The range of T zz over Antarctica is much larger than ±10 E, but here the scale is cut to emphasize various features.(c) The values of the invariant I 1 over Antarctica with complete EIGEN 6C4.(Figure continues on next page).

Figure 11 (
Figure 11 (continued from previous page).(d) The values of the ratio I over Antarctica with complete EIGEN 6C4.Small values of I indicate localities with possible presence of flat bodies (see the theory and the next plot with the strike angle).(e) The values of the strike angle i over Antarctica according to EIGEN 6C4 with the full resolution available (in red its direction to the west, in blue to the east of the meridian).Lambert Glacier (LG), Dehmann Glacier (DG) shown by arrows.(f ) The values of vd over Antarctica according to EIGEN 6C4 with reduced resolution (for more details see figures of vd for the individual segments and places below).(g) The bedrock topography in Antarctica (according to Bedmap2) with a resolution reaching 1×1 km (according to Fretwell et al. [2013]), with the numbers for the names of the main topographic features: 1 -Ross shelf glacier and Ross Island, 2 -Victoria Land, 3 -Transantarctic Mountains, 4 -Wilkes Land, 5 -Marie Byrd Land, 6 -Queen (Dronning) Maud Land, 7 -Princess Elizabeth Land, 8 -Gamburtsev Mountains (under the ice).

Figure 12 .
Figure 12.Segment 1: (a) T zz , (b) I 2 , (c) the RAT14 combination of the gravity from the EIGEN 6C4 and the topography from the Bedmap2 in T zz , (d) the Ross Island topography detail.

Figure 13 (
Figure 13 (continued from previous page).(e) The bedrock topography from Bedmap2, (f ) a zoom of the topography from Bedmap2 for Mt.Sidley, (g) Mt.Sidley -a detail by © Google Earth.

Figure 16 (
Figure 16 (same column, top to bottom).Segment (in a zoom of Figure 15a): (a) the topography from Bedmap2, view with numbers showing the positions of the following details (further and much more detailed zooms) for (b)-(d), 1, 5 and 6, all on land under the ice.The heights in (a) are in meters above sea level.

Figure 17 (
Figure 17 (same column, top to bottom).Topography according to Bedmap2 for the area of lake Vostok LV, (a) on left, and (b) T zz , on right.Note a possible volcano at the centre of the figures (encircled) and the correlation between (a) and (b).