Geophysical Challenges for Future Satellite Gravity Missions: Assessing the Impact of MOCASS Mission

The GRACE/GRACE-FO satellites have observed large scale mass changes, contributing to the mass budget calculation of the hydro-and cryosphere. The scale of the observable mass changes must be in the order of 300 km or bigger to be resolved. Smaller scale glaciers and hydrologic basins significantly contribute to the closure of the water mass balance, but are not detected with the present spatial resolution of the satellite. The challenge of future satellite gravity missions is to fill this gap, providing higher temporal and spatial resolution. We assess the impact of a geodetic satellite mission carrying on board a cold atom interferometric gradiometer (MOCASS: Mass Observation with Cold Atom Sensors in Space) on the resolution of simulated geophysical phenomena, considering mass changes in the hydrosphere and cryosphere. Moreover, we consider mass redistributions due to seamounts and tectonic movements, belonging to the solid earth processes. The MOCASS type satellite is able to recover 50% smaller deglaciation rates over a mountain range as the High Mountains of Asia compared to GRACE, and to detect the mass of 60% of the cumulative number of glaciers, an improvement respect to GRACE which detects less than 20% in the same area. For seamounts a significantly smaller mass eruption could be detected with respect to GRACE, reaching a level of mass detection of a submarine basalt eruption of 1.6 109 m3. This mass corresponds to the eruption of Mount Saint Helens. The simulations demonstrate that a MOCASS type mission would significantly improve the resolution of mass changes respect to existing geodetic satellite missions.


Introduction
Variations of the gravity acceleration on the Earth surface occur on a wide spectrum of time and spatial scales. The complementary gravity missions GOCE (Floberghagen et al., 2011) and GRACE and GRACE Follow On (GRACE/FO) (Tapley et al., 2019), successfully mapped the static gravity field and its time variations, giving a data set for delineating geologic structures of interest and mapping subsurface mass variations globally. The static field was obtained from GOCE with high spatial resolution and the dynamic field from GRACE/FO has lower spatial resolution but high sensitivity.
For instance, water movements in wide inland basins, as the Amazonas (Han et al., 2010), as well as over the oceans (Jeon et al., 2018) and glaciers (Zemp et al., 2019) were successfully detected and studied from the time variation of GRACE data: the analysis revealed the presence of both periodic annual components and long term gravity trends connected to the underground water flow (Döll et al., 2014). GRACE and GOCE satellites were also able to observe the gravity field variations due to the co-seismic and post-seismic deformations of large earthquakes, such as the Sumatra-Andaman (Han, 2006) earthquake or the Tohoku event (Cambiotti, 2020;Han et al., 2015). Plurennial trends in the gravity field variations have been recorded and studied in the Himalaya region (Braitenberg & Shum, 2017;Matsuo & Heki, 2010;Sun et al., 2009) and interpreted as due to the combination of hydrologic, isostatic and tectonic effects.
Cross continent correlation of the geologic structures is now possible (Braitenberg, 2015;Köther et al., 2012) due to the homogenous coverage of observations and the overcoming of geodetic issues that affect the terrestrial measurements. Typically, the problems were connected to the difference in the geodetic reference systems between countries that caused bias at long wavelengths.
The satellite observations from GRACE of the gravity changes have been validated with terrestrial superconducting gravimeters and have shown the complimentarity between terrestrial and satellite observation from GRACE Weise et al., 2009Weise et al., , 2012.
These past satellite missions employed a combination of electrostatic accelerometers, which recorded both inertial and gravitational accelerations; the GNSS tracking of the satellite allowed removing the inertial components unveiling the pure gravity attractions. Recently, the advent of Cold Atom Interferometry (CAI) (Carraz et al., 2014;Ménoret et al., 2018), which offers better noise performance with respect to the electrostatic accelerometers, opened a new era for the satellite observations pushing forward the possibilities to monitor geophysical phenomena at smaller scales.
A feasibility study of a new satellite mission concept utilizing a CAI sensor, instead of the classic electrostatic accelerometers, was the subject of a project that involved Politecnico di Milano, University of Firenze, University of Trieste and that was funded by the Italian Space Agency (ASI).
In the frame of this project, named MOCASS, the three groups analyzed three different aspects of the mission: the University of Firenze studied and defined the characteristics of the instrument for a satellite mission; the Politecnico di Milano performed the mission simulations producing spectral error curves while the University of Trieste group assessed the possible scientific impact.
An overview of the group activities and the overall results has been already published by Migliaccio et al. (2019). In the present work we aim to give more details on the geophysical impact of the mission, presenting in more detail the simulations conducted and some more examples not discussed by Migliaccio et al. (2019). A companion paper to this (Reguzzoni et al., 2021), details the geodetic aspects presenting the simulations of the mission.
We are interested in assessing the magnitude of the gravity signals induced by the redistribution of mass due to hydrological, glaciological and tectonic effects.
We will focus on the High Mountains of Asia area in particular, which is the world largest Mid Latitude glacierized region and is experiencing an important deglaciation (Gardner et al., 2013). In addition to this, the area is characterized by an intense tectonic activity, along the Main Himalayan Thrust, that is responsible of an uplift of the frontal chain and to significant hydrologic variations.

Data and Methods
The area selected for the simulations is in Central Asia and comprises the High Mountains of Asia (HMA) region, including the Himalaya and the Tibet plateau, and part of the Indian Ocean. The HMA region constitutes an important system of the earth where several climatic and geodynamic phenomena take place and influence each other. Such phenomena concur in generating huge mass movements, which involve different spatial and temporal scales. In our simulations we will focus on glacial mass variations, crustal deformation due to tectonic processes, hydrology and volcanic eruptions associated to seamounts. All these phenomena generate mass changes and hence gravity signals over a wide spectrum of temporal and spatial scales: from daily-weekly temporal scales for the volcanic activity to long period trends in which the tectonic and climatic signals superpose.
For estimating the gravity effects, some constraints on the spatial distribution of the volume and density temporal variations are required. We derived most of them from publicly available databases which are presented in detail in the following.

Glaciers Mass Variations
For what concerns the estimate of the glacial effect we focus on the HMA system (Gardner et al., 2013), which comprises the Himalaya, the Tibet and the northern part of the Himalaya chain (Tien Shan). This orogeny is the greatest worldwide in terms of glaciers surface (118,200 km 2 ) and is currently losing glacial masses at a rate of 220 kg/m 2 /yr (Gardner 2224 T. Pivetta et al. Pure Appl. Geophys. et al., 2013). Glaciers in the region have been studied from combination of satellite products including satellite imagery and previous gravity mission as GRACE. GRACE offered constraints on mass fluxes and these data have been used to provide estimates of the deglaciation as testified by the studies of Matsuo and Heki (2010) and Jacob et al. (2012). In order to construct realistic synthetic models and simulate the gravity change we need constraints regarding: (1) The outline of the glaciers (2) The yearly rate of change (dh/dt) of the glacier thickness (3) A DEM (Digital Elevation Model) model to extract the height information of where the glaciers are located For the outline of the glaciers we rely on the Randolph Glacier Inventory (RGI) (Pfeffer et al., 2014) which is a worldwide database of the glaciers. The database is continuously updated and it is distributed in a shape file that contains the polygons of the outlines and the relative attribute table. The polygons have been extracted mainly from the satellite missions such as Landsat or ASTER (Pfeffer et al., 2014) through automatic algorithms. For many glaciers there are more than one outline, which refer to different times of analysis; using a GIS software we extracted the polygons pertaining on a specific time. From the Randolph catalog we extracted the outlines of the glaciers for the whole Himalaya-Tibet area, following the macro-areas outlined in red in Fig. 1a.
The yearly rate of change (dh/dt) of the glacier thickness is derived from the study of Gardner et al. (2013), which evaluated the contribution of the deglaciation to the sea level rise, and provided a detailed analysis of the worldwide ice masses variations. In this study the authors gave an estimation of the yearly elevation change (dh/dt) of the glaciers in the Himalaya region, by analyzing the data from IceSat for the period 2003-2009. The authors reported the average dh/dt values for the macroregions of the RGI catalog shown in Fig. 1a. In the same figure the glaciers outlines are colorized with a color scale proportional to the deglaciation rate of Gardner et al. (2013). Blue values show high deglaciation rates dh/dt [ 0.5 m/yr, yellowish colors show stable glacial thicknesses, while red values, which are found only in one sector of the western Tibet, report glacial thickening.
The last constraint is the DEM model, which is necessary to correctly locate in space the glacial masses. We take advantage of the ETOPO1 model (Amante & Eakins, 2009) which offers global topography/bathymetric values at a resolution of 1'.

Tectonic Deformation and Mass Movements
The tectonic forces cause horizontal and vertical movements of the lithospheric plates that are captured by the GNSS network data. For example in the Himalayan region the convergence of the Indian and Asian plates, that causes frequent earthquakes, is accommodated by the mountain rising which peaks up to 5 mm/yr (Fu & Freymueller, 2012;Liang et al., 2013). This deformation causes mass redistribution and hence a variation in the gravity field. However, such deformation does not involve only the most superficial crustal layers but it also affects deeper crustal units, probably down to the Moho discontinuity.
Two end member mechanisms could be hypothesized (Braitenberg & Shum, 2017;Chen et al., 2018;Yi et al., 2016) as a crustal response to an uplift of the topography: (1) the Moho is shifted upwards by the same amount of the topography and the empty space at the base of the crust is filled by the mantle upwelling. This mechanism is equivalent to the crustal response to unloading (Fig. 1d) (2) the Moho thickens as an isostatic response to the increased topographic load (Fig. 1e) In order to evaluate the gravity effects of these two models we need as constraints the yearly vertical rates due to the tectonic deformation, and additionally a DEM model. The vertical rates are derived from two published studies (Fu & Freymueller, 2012;Liang et al., 2013). Fu and Freymueller (2012) analyzed the GNSS time series data of over 30 stations in the Nepalese area while Liang et al. (2013) used over 750 GNSS stations to represent the current crustal motion and deformation for a wider area Vol. 178, (2021) Geophysical Challenges for Future Satellite Gravity Missions including the Tibetan Plateau, Yun Gun and Sichuan regions. Both the studies analyzed continuous GNSS data (Liang et al., 2013 also included some GNSS stations acquired in Campaign mode) and removed the seasonal periodic oscillations obtaining finally the vertical rates adjusting a linear trend. For the Himalaya area, the comparison of the rates along 2D profiles shows that the observed values are compatible in terms of magnitude with those predicted from a simple model of inter-seismic strain accumulation. As we see in Fig where the uplift easily exceeds 5 mm/yr; the Tibet is generally uplifting as well but at lower rates and with some localized subsidence zones in correspondence of small basins. The zones in the Eastern part of the plateau, are characterized by marked uplift in the northern section, while in the south, towards the Yun Gui plateau, stable terranes or subjected to subsidence are observed.

Large Scale Superficial Hydrology
The hydrologic effect of the Himalaya region was estimated using the GLDAS (Global Land Data Assimilation Systems) products. GLDAS is a project that offers several databases of the land surface state and fluxes in the whole world, at different spatial and temporal resolutions (Rodell et al., 2004). Surface temperature, snow cover, soil moisture, vegetation index are examples of the physical parameters that could be extracted. Each of these is provided at monthly, daily and also 3 h temporal resolutions, while 0.25°and 1°are the spatial resolutions available. The database is constructed incorporating direct observations of the various phenomena from satellite or from ground based observations and modelling procedures. The moisture of the first 2 m of the ground is one of the available physical parameters offered and the most interesting for our purposes.
The moisture database from GLDAS has been already used for correcting geodetic time-series from hydrologic effects: some applications with GRACE could be found in Matsuo and Heki (2010) or for GNSS in Yi et al. (2016). 3 h resolution models are also largely employed for correcting Superconducting Gravimeter observations (Mikolaj et al., 2016). We exploited the GLDAS soil moisture product for the India-Tibet area, with a spatial resolution of 0.25°a nd temporal resolution of 1 month. We selected a period for the analysis that spans from 2002 to 2016. The moisture is provided in various NetCDF grids in terms of equivalent water mass in kg per m 2 .

Seamount Growth
Sea mounts are submarine mountains that are very narrow and very tall, emerging for several thousand meters above the bathymetry. They are formed through volcanism and some are active volcanoes of the ocean floor. A static sea mount bears hazard for shipping when its top reaches shallow depths of about 100 m. The international database of seamounts relies greatly on the inversion of the gravity field recovered from satellite altimetry. Models of temporal changes in the gravity field from satellite altimetry are unavailable, and this technique seems to be difficult to be used for generating gravity changes in time. The complete coverage of the gravity changes is not possible with an altimetric satellite placed on a repeat orbit, since the areas between the orbit paths are not monitored. An alternative could be the geodetic gravity satellites, designed for regular observation of the gravity field, although they have a coarser space resolution. The hazard of a seamount is when it grows to a height that could be intercepted by a ship, taking into account the tallest waves that could affect the area during the strongest wind force and ocean state. In this case the relevance of the monitoring is given by the depth distribution of the top of seamounts below the ocean surface, and the resolution of the mass change that adds mass to the top during a growing episode of the seamount. Multispectral satellite images can be used to detect the seamount once it has reached the surface, but the images are not useful in detecting a growing mount before it reaches the surface. We therefore analyze b Figure 1 a Outlines of the glaciers in the HMA region according to RGI catalogue (Pfeffer et al., 2014), colorized proportionally to the (de)glaciation rate in m/yr. The sub-regions defined in the RGI catalogue are outlined in red and labeled with the following acronyms: HA-P: Hissar Alay and Pamir; TS: Tien Shan; QS: Qilian Shan; S-ET: South and East Tibet; HS: Hengudan Shan; IT: Inner Tibet; W-C-E Himalaya: West-Central-East Himalaya; HK-K: Hindu Kush and Karakoram. In the background the ETOPO1 hill-shade map. b GNSS observations available in the study area; color scale gives the vertical rates. For black dots only horizontal rates are available and do not enter this study. In the background the ETOPO1 hill-shade map. c topo/bathymetry of the area from ETOPO1 with the seamounts location according to the Wessel et al. the relevance of catching the signal of a growing sea mount, next to the sensitivity of documenting existing ones.
The sensitivity analysis of the gravity signal due to sea mounts is done in the area of the Indian Ocean. Our estimate of the signal of a growing sea mount is not tied to the specific area of the Indian Ocean, but is intended to be an estimate in general terms. Therefore, the question whether specifically in the Indian Ocean the seamounts are volcanically active is beside the point, since we analyze this area due to the large amount of seamounts, generating a reliable statistic. Figure 1c shows with black squares the seamounts distribution according to the database of seamounts compiled by Wessel et al. (2010) in which the seamounts were identified using the ERS-1/Geosat gridded vertical gravity gradient. The database contains longitude and latitude of each seamount, radius and height, and the age of underlying seafloor. Figure 1f shows the statistical distribution of the depths of the tops of the seamounts; we derived the tops depth from the ETOPO1 mode. As it is evident most seamounts have shallow tops, mostly below 100 m.
For our scopes we are interested in determining the depth of the top of the seamount with the aim of simulating the eruptions, just imposing different extra volcanic masses on top of it. The volumes of the eruptions tested are derived from estimates of historical eruptions, as we will discuss in the following, and they are in the range of 10 8 -10 9 m 3 .

Forward Modeling of the Gravity Fields
For calculating the temporal gravity variations we discretize the masses of the various geophysical phenomena through tesseroids (Uieda et al., 2016) except for the seamounts, in which we exploit a numerical solution of the Poisson equation through a finite element solver. The tesseroids are the equivalent of prisms in a spherical reference system and like prisms represent a 'building block' through which to build arbitrary complex mass distributions. In order to calculate the gravity effects of glacial thickness variations we select the outlines from the RGI catalogue pertaining to the HMA region and for each outline we extract the topography quotas from the ETOPO1 model (Amante & Eakins, 2009). For each topographic point we construct a tesseroid with a basal area of 0.0167°9 0.0167°, matching the ETOPO1 resolution, and a thickness corresponding to the yearly height change of the glacier, according to Gardner et al. (2013). The density of each tesseroid is fixed to be equal to 900 kg/m 3 . The yearly mass change per unit area is reported in the supplementary Figure S1a.
The gravity effect of the glaciers for the Himalaya region is calculated on a regular grid of 0.25°resolution at 250 km of quota and it is shown in Fig. 2a.
For the tectonic signals we proceed similarly as for the glacier signal; first we interpolate the sparse GNSS observations over a regular grid of 0.25°resolution for the whole area. We then calculate the gravity contribution of the ETOPO1 topography subjected to the interpolated vertical movements, assuming the crustal density q c to be constant and equal to 2670 kg/m 3 . The mass model is reported in Figure S1b. The crustal response is then assessed according to the models presented previously: in case of a pure crustal uplift the Moho is shifted up by the same amount of the topography and it is replaced by the denser mantle (q m ¼3200 kg/m 3 ). Conversely when subsidence occurs the Moho is shifted downwards, and light crust now replaces the mantle. In both the situations the signals due to the vertical movements occurring on surface are amplified by the deep crust mass redistribution. In our case the Moho interface in the HMA region is approximated by an isostatic Moho. The signals for the crustal uplift model are shown in Fig. 2b.
The isostatic response for a surface uplift predicts a lowering of the Moho interface, with sequent increase of the crustal thickness, replacing the dense mantle. Here the superficial mass excess is counterbalanced by a mass loss at the Moho level. Given the densities of the crust and mantle we could calculate the crustal thickness variation Dr for a surface mass variation on topography (q c Dh) by the following equation: The gravity variations associated to the isostatic effect are shown in Fig. 2c. We see they are in 2228 T. Pivetta et al. Pure Appl. Geophys. general one order of magnitude lower than the pure crustal uplift response. The last effect we inspect for the HMA is due to the surficial hydrology which is estimated relying on the monthly GLDAS soil moisture product. Each monthly grid is discretized in tesseroids with resolution equal to that of the model (0.25°) and the gravity effect is calculated at 250 km obtaining a volume of data g z (lat, lon, time). For each longitude and latitude couple we extract the g z time series and we fit it with the following equation: obtaining a map of the linear trend (m) and of the yearly hydrologic oscillation amplitude defined as ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . Both the long period and seasonal mass changes are reported in Figure S1c and d. Figure 2 Simulated gravity fields. a yearly gravity change due to glaciers thickness variations. b tectonic yearly gravity change due to crustal uplift. c yearly gravity change due to an isostatic response of the crust to surface vertical movements. d hydrologic long period trend from GLDAS. e annual signal amplitude due to hydrology. In all figures the purple outline reports the 3500 m topographic isoline and black line shows the coastline; the various thick lines that cross the anomalies are the traces of the profiles used for the spectral analysis discussed in the following sections. f gravity change due to a seamount growth as a function of the maximal thickness of the eruptions. g cross section of the seamount geometry Vol. 178, (2021) Geophysical Challenges for Future Satellite Gravity Missions Figure 2e shows the periodic annual component in the HMA region; the most prominent signal comes from the Ganges and Brahamputra region, where the yearly signal reaches the 0.004 mGal amplitude. Another interesting periodic signal is observable in the NW region of the map, which corresponds to the Fergana Valley (Uzbekistan).
The last Fig. 2f, g show the simulations of a sea mount eruption and the resulting gravity change. Mechanically the eruption should be a non-explosive magma outflow of basaltic lava at the top, which is immediately cooled by the ocean, analogous to the formation of pillow lavas. We consider a seamount with a cone-like shape of radius 12,500 m, height 5800 m, top at 200 m below sea surface, as shown in sketch 2 g. We successively add magmatic volume to the top, with an increasing height of 5 m, 10 m, 50 m, 100 m, 120 m, 140 m, 150 m and 160 m. Density of the seamount is 2800 kg/m 3 . For each magmatic event the gravity field change is displayed in Fig. 2f with color code proportional to the maximal thickness of the eruption.
To further prove the goodness of our simulations we compare the summed effects of the tectonic and glacier signals with the yearly gravity trends observed by GRACE. We rely on the Jet Propulsion Laboratory (JPL) monthly solutions synthetised for the period 2002-2016, corrected for the hydrologic component estimated through the GLDAS products. From the various monthly grids we extract the long period trend which is shown in Fig. 3a. In the simulated fields (Fig. 3b) the high frequency contribution is filtered out through a low-pass Gaussian filter, with cut off wavelength of 500 km.
From the comparison of the two plots we note a general good agreement both for what concerns the pattern of the anomalies and their amplitudes. For instance, the arch shaped anomaly due to the glaciers is present in both observation and models in the southern part of Tibet and in the Central Himalaya. In the western sector (Lon: 75°; Lat: 30°) GRACE predicts a wider anomaly zone compared to our simulations. This is probably due to a long-term water mass variation associated to intense anthropogenic water extraction (Saji et al., 2020) which is not modelled in GLDAS products. The crustal tectonic uplift seems to be able to explain part of the positive trend in the Central-Eastern Tibet, the mismatch in amplitude could be caused by other components not modelled such as water volume variations associated to the endorheic lakes, as suggested by Zhang et al. (2017).

Geophysical Impact Assessment
The final step of the research aims at setting the results of the geophysical simulations of the expected gravity field into the context of the error curves of the innovative satellite gravity mission (Migliaccio et al., 2019;Reguzzoni et al., 2021). In order to evaluate the significance of the satellite mission using the atom gradiometer on board, we estimate the percentage improvement of the completeness of phenomena that can be observed. Given a topic of investigation, as for instance detection of the mass loss of a glacier, it is obvious that an increased precision and accuracy of the gravity observation and a higher spatial resolution, allows to detect smaller glaciers loosing minor ice mass due to melting. The question we pose is the significance of the improvement, which is based on the quality jump of the investigations that can be achieved only with the new mission. One criterion starts from the statistical evaluation of the numerosity-size distribution of the glaciers, or the other phenomena we have investigated (hydrologic basins, sea mounts, mountain dynamics). Given a detection threshold, the percentage of the glaciers that can be monitored is a criterion that we estimate. In the following the results for mountain building process, hydrologic basins, glaciers and sea mounts are discussed.

Time Variable Phenomena
First, the gravity change of the sum of the glaciers is analyzed. With the scope of comparison with the spectral error curves, it is necessary to determine a gravity change value corresponding to the characteristic wavelength. This is done by tracing a series of orthogonal profiles across the mountain belts, along which the gravity change rate variation is determined. We proceed to estimate the typical wavelength and amplitude of the considered phenomena. For each profile reported in the previous Fig. 2 we fitted the following function according to Jekeli (1981): The function resembles a Gaussian function assuming a [0 and for 0\x \ p. Jekeli (1981) reported that a is related to the inverse of the second moment of the Gaussian distribution (r 2 ). We estimated the wavelength of the signal to be 4r; b is the amplitude of the signal. The following Fig. 4a reports an example of the fitting procedure over a profile traced on the glacier simulation and the fitted curve from Eq. 3 (fitting curve in red).
The estimated wavelengths and amplitudes for each profile are reported through dots in the error curves of the mission. A phenomenon can be detected if the signal (dot) lies above the noise curve.
In the specific case of the glaciers in Himalaya and the other mountain ranges surrounding the Tibetan plateau, the glaciers are in close proximity, and coalesce to form a broad glacier surface, producing a long wavelength signal. It is a more favorable situation than other areas with single glaciers, as could be the case of the Alps, with a much lower density of glaciers. The wavelength of the signal of a single glacier, approximately a point mass, is large at the height of the satellite, for instance at 250 km, and is in the order of degree 50 in the spherical harmonic expansion. As seen in Fig. 4b the present deglaciation is well above the error curves of both MOCASS and GRACE. Of interest is the observation of variations in the average yearly rate, for instance to determine whether recent variations of accelerated deglaciation are visible. Here MOCASS shows to be superior, because it is able to catch lower fluctuation rates (less than 50%; green dots), which instead lie on the error curve of GRACE, and therefore cannot be detected.
Next we consider the gravity signal at 250 km height due to the changing mass of grouped glaciers of increasing areal extent, imposing three different deglaciation rates (-0.12, -0.40, -0.89 m/yr), which are realistic values derived from Gardner et al. (2013). The characteristic wavelength is found through fitting a Gaussian curve, of which the half width (sigma) is set equal to a fourth of the signal wavelength. The amplitude of the Gaussian curve is taken as the amplitude of the estimated signal at the specific wavelength. We compare this gravity change value with the spectral error curves of three satellites, the GRACE, GOCE and MOCASS in two configurations, estimated for one year of data acquisition. The two MOCASS error curves are for the MOCASS nadir-pointing double-arm configuration, with atomic single-arm or double-arm gradiometer (Migliaccio et al., 2019), at low orbit of 239 km height. The comparison assumes that the mass change of glaciers is steady state, so a yearly sampling of the signal is sufficient. The results are shown in plots 5b and 5c. Due to the calculation height of 250 km the wavelength of the signal, also for a small glacier, is relatively large and is in the range of degree 50 and Vol. 178, (2021) Geophysical Challenges for Future Satellite Gravity Missions 2231 70. This because the gravity signal also for a point mass widens with increasing height. The color codes of the dots correspond to four glacier areas, in which smaller glaciers have been gathered to form a compound glacier of increasing size. The four glacier areas are shown in Fig. 5a. The glacier areas smaller than 500 km 2 are below resolution of the satellites.
Starting with an area of 1900 km 2 the deglaciation can be observed, for rates above 0.40 m/yr, but only with MOCASS, not with GRACE. Glaciers as big as 4750 km 2 must have a rate of at least 0.4 m/yr to be seen by GRACE, while also smaller rates are detected by MOCASS. A further question is, whether the deglaciation signal of two neighboring glaciers of the HMA can be resolved, and which is the smallest resolvable distance between the glaciers. We estimate this distance by approximating two glaciers each by a point mass. The volume change is equal to the area of the glacier times the assumed deglaciation rate of -0.40 m/yr. Given the spectral error value dg(n) of the satellite, two glaciers can be resolved if the difference between the gravity of one glacier and the gravity minimum value of the superposition of two glaciers is smaller than the degree error (Fig. 5d). The distance between the two peaks of the point masses is the lateral resolution. The distance dist is given by the following equation, with G the gravitational constant, M the yearly mass loss of the two glaciers,h the calculation height, and r the error at degree n.
Given the calculation height and spectral error, and making the area of the glacier vary between 500 and 12,000 km 2 , we determine the smallest distance between the two glaciers that can be resolved. The deglaciation rate is fixed to 0.40 m/yr (Fig. 5e). The smaller glaciers can only be resolved at greater distances. At degree n = 50, corresponding to a distance of 792 km, MOCASS resolves two glaciers of about 800 km 2 size, whereas GRACE can resolve two glaciers of more than double size, of about 3000 km 2 . At degree n = 100, MOCASS resolves glaciers of size 900 km 2 , GRACE only glaciers of size 12,000 km 2 .
Then, we compare the signal amplitudes with the error curves for the other phenomena considered. In all cases the profiles, for which the traces are reported Figure 5 a Four classes of glacier areas of increasing areal extent in which glaciers have been aggregated. Outlines mark classes: blue: 100 km 2 , green: 500 km 2 , red: 1900 km 2 , and the total area of the glaciers in the window is 4750 km 2 . In purple the glaciers outlines from RGI. b Spectral comparison of the simulated gravity change for single glaciers or groups of glaciers in HMA with the error curves of satellites GRACE, GOCE and MOCASS in two configurations. The error curves correspond to one year of gravity observation. c area of glaciers and the gravity change in one year due to height loss of the glaciers. Different rates of height loss are assumed (-0.12 m/yr, -0.40 m/yr, -0.89 m/yr). The single simulations are compared to the error curves in the plot b of the same figure. The gray and yellow areas report the deglaciation rates and glacierized areas sensed by GRACE and MOCASS, respectively. d Minimum distance between two gravity peaks to be distinctly detected. in Fig. 2, are fitted with a Gaussian function to find the characteristic wavelength and its amplitude, for comparison with the spectral error curves. We start with the tectonic signals, which are shown in Fig. 6a. The greater part of the signal can be resolved both with GRACE and MOCASS, but some features are completely lost in the GRACE observation, as the dot at n = 80 which corresponds to a localized subsidence in the North Tibet plateau. To understand the origin of the subsidence, the knowledge of the gravity change rate is necessary, which can be resolved by MOCASS.
The hydrologic signals sensitivity is shown in the plots 6b and 6d. In both cases, aperiodic and periodic Figure 6 a Yearly gravity change rate at 250 km calculation height for Tibet-Himalaya due to tectonic uplift measured by GPS. The black dots display the characteristic wavelength and amplitude for the profiles in Fig. 2b. The satellite error curves are 2-years solutions. b Hydrologic linear trend in the HMA region at 250 km calculation height. The dots display the characteristic wavelength and amplitude of the signal along the profiles shown in Fig. 2d. Green dots are for Pakistan area; red dots are relative to Eastern Tibet. The satellite error curves are 1 year solutions. c Gravity change at 250 km calculation height for the magmatic growth of a seamount compared to the satellite error curves. The labels report the height increase of the submerged volcanic edifice. Monthly error curves for different satellite missions are shown. d Hydrologic annual oscillation amplitude in the HMA region at 250 km calculation height. The dots display the characteristic wavelength and amplitude for the profiles shown in Fig. 2e patterns, two areas are analyzed by fitting a Gaussian curve to the anomaly profiles. For the trends (Fig. 6b) we observe that the signals are very close to the error curves; GRACE performs better at low spatial frequencies (low n) while MOCASS, with the flat error curve for n [ 40, detects some high frequency signals that are below the GRACE error curve. This simulation shows the complementarity of the missions; a combined gravity field model including GRACE follow-on and MOCASS could greatly improve the understanding of the dynamics also in small basins. The signal generated by submarine mass variations due to seamount growth sensitivity is plotted in Fig. 6c. Major eruptions can be seen both in GOCE and GRACE, but only MOCASS is sensitive to the more realistic eruptions, increasing height by a few meters. In the case we chose, a height increase of 5 m would be seen by MOCASS, whereas GRACE would need an increase of more than 10 m. In terms of volumes, for an eruption of 3 m, the erupted volume is 4.9 10 8 m 3 , of 5 m, the erupted volume is 8.2 10 8 m 3 , and for 10 m, the erupted volume is 1.6 10 9 m 3 . These numbers can be compared with the volume erupted with the Saint Helens volcanic eruption, which produced 1.6 10 9 m 3 of volume in a single eruption. Our estimate is of first order and could be refined in more detail by including a model of the magmatic production of a sea mount, with the feeder channel from lower crust and mantle and the magma chamber. In the case that the magma system outpours the magma, but is not refilled from below, the mass change would be smaller than we have estimated. Our estimate assumes that the erupted volume is replaced in the magma chamber through rising magma in the feeder channel.
Finally, we compare the amplitude and wavelength of the annual signal, calculated over different profiles (traces are reported in Fig. 2e), with monthly satellite error curves in Fig. 6d. The amplitudes and wavelengths of the phenomena have been calculated with the approach already presented in the previous sections.
The signals are all well above the noise threshold of both GRACE and MOCASS missions.

Glaciers Static Detection
The final analysis on glaciers regards the problem of total mass detection of a glacier. This is relevant when the area of the glacier has been determined from remote sensing as with SAR and multispectral images, but the thickness is unknown (e.g. for Julian Alps Colucci and Ž ebre (2016)). In the simulation we calculate the gravity field generated by the ice-mass. It may be argued, that if the glacier surface is known, and the topography erroneously has been set equal to this glacier surface, the glacier mass must be calculated using the density contrast of ice relative to the rock constituting the topography. To make the analysis more general, and assuming the topography represents the bottom surface of the glacier, we use the density contrast of ice against air. In the case that the contrast must be calculated as ice against rock, the gravity signal will be greater by a factor 1.6. This because the density of ice is close to 1000 kg/m 3 , the density contrast of ice against rock is -1600 kg/m 3 , assuming a granite rock density of 2600 kg/m 3 . We calculate the glacier volume starting from the areal extent using an empirical formula for the glacier height (Bahr et al., 2015), with S the area in km 2 and V the volume in km 3 : Single glaciers of different size are selected for the HMA from the Randolph database, and discretized into tesseroids of 1' sides. The base of the tesseroids is set equal to the ETOPO1 topographic height, the top is calculated through the empirical height formula. Gravity calculation height is 10 km, so as to be above the topography. Only for the smallest area of 10 km 2 a single tesseroid of square footprint is used. Again the signal is fit with a Gaussian curve to determine wavelength and amplitude, and these values are compared to the 5 years error curves of the satellites. Approximately the observation thresholds are a glacier area of 250 km 2 for GRACE, 100 km 2 for GOCE, and 10 km 2 for MOCASS, as seen in Fig. 7.
Considering the statistical distribution of the number of glaciers versus size, MOCASS achieves a significant step, because at sizes around 100 km 2 the number of smaller glaciers jumps significantly Vol. 178, (2021) Geophysical Challenges for Future Satellite Gravity Missions with respect to larger glaciers from 20 to 80 glaciers, as seen in Fig. 7. With smaller areas, close to 10 km 2 the number of glaciers is 20 fold, which demonstrates that MOCASS achieves a very significant increase in glacier monitoring, because it crosses the border between copious and scarce number of glaciers, at least for the HMA. MOCASS covers therefore a much larger total mass of glaciers, including the small glaciers, against GOCE and GRACE, which catch the larger glaciers, but catch the signal of small glaciers only as aggregate signal, when the glaciers are close to each other and coalesce to one signal. Sparse and small glaciers are not observable by GOCE and GRACE.

Discussion and Conclusions
First estimates show that MOCASS brings significant improvement in the gravity field with respect to GRACE. An operational use of a satellite mission for glacier monitoring in mid-latitudes can presently not be made with GRACE because the resolution is too coarse. To catch 60% of the glacial areas in HMA a resolution of 10 km 2 is needed, provided by MOCASS. The signal would provide the mass loss of the integrated glaciated area, as well as the possibility to define present ice volumes starting from the glacier area detected from Sentinel 1 and 2, either through SAR or multispectral imaging.
Tectonic vertical movement is caught by GRACE, but the local anomalies, as in the Yunnan province, can only be caught with MOCASS. With gravity observations the geodynamic reasons for the subsidence can be distinguished, as tectonic stresses and sediment compaction.
For studying the hydrology of Tibet, only MOCASS can detect the predicted changes from the GLDAS, as the observed reduction in moisture found in Eastern Tibet. In Eastern Tibet GRACE has a positive gravity increase rate, that does not match the moisture decrease. Here the observed GRACE signal b Figure 7 Gravity signal generated by the mass of different-sized glaciers in HMA. a the location of the selected glaciers and their outline. b gravity signal compared to error curves. The glacier height is estimated from the areal extension through an empirical formula. The signals are filtered with Gaussian low pass filter (500 km cutoff wavelength). Only the small glacier with 10 km 2 is also filtered with 250 km. The error curves are for 5-year observation time. c Statistical distribution of number of glaciers in HMA with respect to their size. The horizontal lines show the detectability ranges of the three satellites GRACE, GOCE and MOCASS. Blue: histogram of glaciers in function of glacier area; red: Complementary cumulative distribution function (ccdf). Y-axis: left for histogram, right for ccdf. The first bin of the histogram reaches the value of 1900, the second about 700 2236 T. Pivetta et al. Pure Appl. Geophys. probably contains the superposition of the tectonic and hydrologic mass change of opposite sign, and the two effects cannot be resolved. A mission as MOCASS would allow to separate the two hydrologic and tectonic effects in this strategic area.
In the Alps, where the hydrologic basins are much smaller (4.5 10 4 km 2 Po Basin vs. 2.5 10 6 km 2 of the Brahmaputra-Indus-Ganga systems), the resolution of GRACE is far from the expected signal, and only MOCASS would allow to detect the hydrologic mass changes.
A further innovation with the MOCASS satellite would be the continuous monitoring of mass changes at the top of sea mounts. The appearance of an uncharted ocean island is a rare event, due to volcanic submarine activity. Although a rare event, in case of ship collision the risk is great, be it ecological or of human loss. For submarine vehicles also the growth of a deep top of a seamount is a danger, whereas for ships a shallow top is hazardous. The seamount with a shallow top below the sea surface, is not observed from optical observations as obtained from Sentinel 2 or SAR observations obtained from Sentinel 1. It may be detected from satellite altimetry, but the repeat pass orbit has uncharted areas, in which a growing seamount could be fortuitously placed. The MOCASS mission would detect a mass change of a seamount of 12 km radius, adding a volcanic mass to the top of 10 m, which would result in a mass equal to the erupted volume of Mount Saint Helens.
The above examples have been tested in their first orders, and more study is required to underpin the results before deciding on such a mission. One aspect to be studied in greater depth is the method of comparing a local effect, as the one we have synthesized, with the global error curves derived from spherical harmonic expansion. We propose to make a thorough study, both transforming the local fields in the spectral space, as also comparing the local synthesized fields with examples of noise simulations of the MOCASS satellite.
The further studies that could be made are an in depth study of the four geophysical phenomena we have considered up to now, and the study of the added value of a gravity mission to detect the different cycles of a seismogenic fault. In the following, we mention the further study that we think is relevant for the phenomena discussed in this work.
Glaciosphere: investigate with which methods the glacier signal can be separated in the MOCASS synthetics from the remainder gravity signal. Further, how the volume of a glacier can be effectively determined given the digital terrain model which generally will give the top of the glacier surface, and the outline of the glacier from Sentinel 1 and 2. Determine the effective ability of MOCASS to be an operational satellite for glacier mass and glacier mass depletion detection in High Mountain of Asia and other high mountain areas worldwide.
Hydrosphere: given the discrepancies between the predicted signal from the soil moisture models and the direct observations of the underground water level from wells as shown in Chen et al. (2018) the hydrologic model should be improved, so as to find a perfect match between GLDAS (moisture model) and the wells observations of underground mass changes.
Sea mounts: we have demonstrated that a mass change on top of the sea mount is detectable, in principle. The volcanic model should be improved, containing the volcanologic process of magma eruption, to be sure that the mass erupted on the top is not just shifted upwards from the volcanic edifice, but that effectively there is a mass change in the vertical column.
Our final conclusion is that the MOCASS mission has great potential for monitoring geophysical phenomena, but that there is some more investigation needed to determine with high confidence the applicability of the observed gravity field as a monitoring device.
Author contributions CB was responsible and supervised the project, providing the main conceptual ideas and following the project-stages. TP performed the coding, database retrieval and simulations over the HMA region, DFB performed the simulations of the seamounts eruption. TP and DFB compared the signals with the error curves of the gravity missions. All the authors discussed the results of the simulations and assessed the impact of the mission.

Funding
Open access funding provided by Università degli Studi di Trieste within the CRUI-CARE Agreement. The MOCASS study has been funded under Italian Space Agency contract N. 2016-9-U-0 Proposal of a satellite mission and sensor concept based on advanced atom interferometry accelerometers for high resolution monitoring of mass variations on and below the Earth surface''.

Data availability
The data used in the study was retrieved on public databases, as described in the manuscript.

Code availability
Calculations have been made with Custom codes.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.