1 Introduction

Successful satellite gravity missions such as GRACE (Tapley et al., 2004) and GOCE (ESA, 1999) have provided plenty of data for monitoring the Earth gravity field and its changes. More recently, the GRACE-FO twin spacecrafts (Kornfeld et al., 2019) were launched on May 22nd, 2018, in order to continue the monitoring of the Earth water transport with the aim of detecting changes in underground water storage, in the amount of water in large lakes and rivers, soil moisture, ice sheets and glaciers, and in sea level caused by the water transfer from continental areas to the ocean (Frappart & Ramillien, 2018).

Ongoing studies such as NGGM are also being carried out for the long-term monitoring of the time-variable gravity field with high temporal and spatial resolution (Haagmans et al., 2020). A deeper knowledge of several geophysical phenomena in many interesting areas can be made possible by further and more accurate satellite-based investigations of the Earth composition. As for the Solid Earth, knowledge improvements are expected for megathrust earthquake modeling, study of mass distribution and transport inside a volcano, mass transport at an increased spatial and temporal resolution, study of Moho discontinuity and upper mantle, and so on. Even more significant outcomes are expected in oceanography, hydrology, and cryosphere science. Just to mention an example, GOCE significantly contributed to the knowledge of the mean dynamic topography at global scale (Knudsen et al., 2019) thanks to a geoid estimation with an accuracy of 1–2 cm up to a maximum spherical harmonic degree 200–220 (ESA, 1999); however, space and time resolutions are still insufficient for small and closed or semi-closed basins, like the Mediterranean Sea, and an improved gravity mission could fill this gap. Besides, ice masses carry a significant gravity signal: time-variable gravity solutions allow to quantify the ice-sheet mass variations and their impact on climate. All these improvements can give great benefits to climate and climate change studies (Pail, 2015). The geophysical impact of future gravity missions with significantly improved performances with respect to GRACE and GOCE has been deeply studied, e.g., in the framework of the Global Geodetic Observing System (GGOS) Working Group for Satellite Missions (see again Pail, 2015).

While electrostatic sensors are approaching their ultimate performance, atomic quantum sensors have shown the potential to drastically increase the performance of inertial measurements (Peters et al., 2001; Sorrentino et al., 2014; Tino et al., 2013). For example, the GOCE gravity gradients showed poor performance in the lower frequency band, where the noise power spectral density (PSD) increases as the inverse of the frequency (Rummel et al., 2011). In order to improve the GOCE performance (of the order of 10 mE/√Hz in the 5–100 MHz band, where 1 E = 10−9 s−2), future Earth gravity missions will require gradiometers with sensitivity of the order of 1 mE/√Hz over a wide spectral range. Indeed, the remarkable stability and accuracy that atom interferometers have reached for inertial measurements can play a crucial role for science and technology. The architecture of a quantum gradiometer, e.g., a gradiometer built on cold atom interferometry for the purpose of gravity field recovery, was proposed in Carraz et al. (2014). This new technology is expected to enable high-sensitivity measurements of all diagonal elements of the gravity gradient tensor, as well as the full spacecraft angular velocity.

The technological development and performance of a quantum accelerometer were investigated in a more detailed way in several studies (Douch et al., 2018; Mottini & Anselmi, 2019; Trimeche et al., 2019) and supported by the European Space Agency (ESA) in the context of the “Study of a Cold Atom Interferometry (CAI) gradiometer sensor and mission concepts” (ESA, 2014).

In the years of 2014–2018, the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) started funding the “Collaborative Research Center 1128: Relativistic Geodesy and Gravimetry with Quantum Sensors (geo-Q)”, consisting of 21 projects. The aim of this program was to explore new frontiers and techniques for estimating the Earth gravity field and monitoring the Earth mass redistribution at a new level of accuracy, providing improvements needed for the GRACE and GRACE-FO data analysis. One of the three main research areas of this program focused on the design and development of new instruments and sensors for gravity missions. Particularly, the study of quantum gradiometer was also funded by ESA. During this research, the cold atom interferometry and the optical clock concepts were studied, the use of atomic clocks to directly observe gravity potential differences by measuring the relativistic redshift between clocks (“relativistic geodesy”) was proposed and the contribution of the quantum sensors to improve the Earth gravity field was evaluated (Müller & Wu, 2019).

Between 2016 and 2019, Thales Alenia Space Italia performed the study of the “Cold Atom Inertial Sensors (CAI): Mission Application” for ESA with the support of Marwan Technology, Politecnico di Milano and Technical University of Delft (Mottini & Anselmi, 2019). The purpose of this study was to investigate the possible implementation of a future generation of gravity missions using gradiometers based on the cold atom interferometry. In the scope of this study, the design of a mission with the cold atom gradiometer on board and the achievable system-level performances were investigated.

Finally, it should be remarked that at the IUGG General Assembly in Montreal in 2019, the International Association of Geodesy (IAG) decided to establish a project on Novel Sensors and Quantum Technology in Geodesy (QuGe); one of the main purposes of this project is to investigate the enhanced prospects for satellite geodesy that quantum technologies will open up (Poutanen & Rózsa, 2020).

In this paper, we will focus on the data analysis of the MOCASS mission proposal (Mass Observation with Cold Atom Sensors in Space), carried out by a team of researchers from Italian universities and research institutions in the period from March 2017 to June 2019 under a contract by the Italian Space Agency (ASI), in the framework of preparatory activities for future missions and payloads of Earth Observation. In Sect. 2 the mission concept will be explained, introducing the cold atom gradiometer noise characteristics and mission scenarios for which data simulations and analysis were performed by applying the so-called space-wise approach. In Sect. 3 the results of the numerical simulations will be described and compared to those of GRACE and GOCE missions. In Sect. 4 a discussion on the obtained results and some conclusions will be presented. Regarding the science applications of the results of this study, geophysical investigations on glaciers mass variations, tectonic deformations, mass movements and superficial hydrology in the Himalayan region and seamount growth in the area of the Indian Ocean can be found in Pivetta et al. (2021).

2 The MOCASS Mission Proposal and Data Analysis

The MOCASS study was carried out by a research team consisting of Politecnico di Milano, University of Trieste and AtomSensors srl (a spin-off of the University of Florence). The goal of this study was to determine:

  • the technological characteristics of the cold atom gradiometer capable of acquiring observations of gravity gradients at satellite altitude;

  • the estimation of the global signal characteristics as a result of the data analysis;

  • the requirements on the signal to be detected in order to identify the geophysical phenomena of interest.

An overview of the proposal and results of the study is presented in a paper by Migliaccio et al. (2019).

2.1 The Cold Atom Gradiometer

Future improvements of the GOCE mission concept can be achieved by going beyond the technology of electrostatic gradiometers and taking advantage of a new generation of quantum sensors. In recent years, a new and highly improved gradiometer concept has been proposed (Carraz et al., 2014), exploiting instruments that have been refined and tested over nearly three decades, based on the principle of Cold Atom Interferometry. Such instruments operate as inertial sensors and have been used for fundamental physics experiments. A gradiometer in one given direction (e.g., the radial one) can be obtained by measuring the acceleration of two clouds of cold atoms separated by a certain distance. At the same time, this instrument allows to measure the rotation angle around the axis perpendicular to the plane of motion of the atom clouds. Since this scheme can be repeated in the three directions, it will be possible to determine all the diagonal elements of the gravity gradient tensor and the angular rate vector.

For this type of instruments, it has been proved that a sensitivity of 3.5 mE/√Hz over a wide spectral range can be reached. An interesting feature is that its noise power spectral density is flat at low frequencies in contrast to the spectrally colored noise of the electrostatic accelerometers (see Fig. 1).

Fig. 1
figure 1

Noise power spectral density for the cold atom gradiometer, as compared to the one of GOCE (Carraz et al., 2014)

Full size image

As a consequence of its spectral characteristics, such an instrument might allow to meet the requirements of a mission dedicated to the observation of both the time-variable gravity field (like GRACE) and the static field (like GOCE). Of course, in the former case a longer lifetime is required, and this would reflect in the mission characteristics.

2.2 Space-Wise Solution Strategy

In the case of the MOCASS study, the input data were the simulated second derivatives of the anomalous gravitational potential T at given points along the orbit; the output data were the spherical harmonic coefficients of the Earth potential, estimated from the cold atom gradiometer observations by applying the space-wise approach.

The main idea behind the space-wise approach (see e.g., Migliaccio et al., 2004a; Pail et al., 2011) is to estimate the spherical harmonic coefficients of the geopotential model by exploiting the spatial correlation of the Earth gravity field. For this purpose, a collocation solution can be devised, modeling the signal covariance as a function of spatial distance and not of time distance, as it is done in the case of the noise covariance. In this way, data which are close in space but far in time can be filtered together, thus overcoming the problems related to the strong time correlation of the observation noise.

Although a unique collocation solution would be theoretically clean and desirable, it is not computationally feasible due to the huge amount of data which usually must be processed for this type of satellite missions (Pail et al., 2011). The dimension of the system to be solved would be in fact equal to the number of input data. For this reason, the space-wise approach is implemented as a multi-step collocation procedure (Reguzzoni & Tselfes, 2009), basically consisting of:

  • Wiener filtering (Albertella et al., 2004; Papoulis, 1984) applied to the data along the orbit in order to reduce the highly time correlated noise of the gradiometer;

  • spatial interpolation of filtered data to obtain values over spherical grids at mean satellite altitude, by applying collocation to local patches of data (Migliaccio et al., 2007);

  • spherical harmonic analysis of gridded data by numerical integration (Colombo, 1981) to retrieve the geopotential coefficients.

In the space-wise strategy as implemented in GOCE data analysis, the procedure was iterated until convergence to recover the signal frequencies cancelled by the Wiener filter along the orbit and to correct the estimated rotation from the gradiometer reference frame to the local orbital reference frame (LORF) (Migliaccio et al., 2004b). In the MOCASS simulations, no iterations were required (see Fig. 2) due to the prevailing white behavior of the noise spectrum of the cold atom gradiometer and to the assumption that the spacecraft (and therefore the gradiometer) is kept aligned with the LORF; the latter assumption would be justified by the use of an enhanced spacecraft attitude control system, like the one devised for NGGM (Bacchetta et al., 2017). The scheme in Fig. 2 without iterations allowed for higher speed in computations and hence the possibility to analyze several case studies.

Fig. 2
figure 2

The space-wise approach basic scheme

Full size image

Note that in the space-wise approach a prior model has to be used in order to reduce the spatial correlation of the signal, which is necessary when applying collocation gridding on local patches of data. In other words, after removing the prior model contribution from the observations, low frequency information in the residual signal is significantly reduced, leading to a signal covariance function with a correlation length comparable with the size of the local patches used for the gridding. For this purpose, the gravity field model derived by SST (Satellite-to-Satellite Tracking) data was employed (Migliaccio et al., 2010), thus avoiding the use of external information.

2.3 Numerical Simulations for the MOCASS Study

The numerical simulations for the MOCASS study were performed considering as reference orbits the GOCE orbits from November 2009 to January 2010 (about 2 months of data along a “high” orbit at an altitude about 259 km) and from February 2013 to April 2013 (same time span of about 2 months, corresponding to a “low” orbit at an altitude about 239 km).

EIGEN_6C4 (GOCE/GRACE/ground data combined model) up to harmonic degree and order 360 was used as the reference model (Förste et al., 2014) in order to simulate the spherical harmonic coefficients as input.

Two different pointing modes of the instrument were considered: nadir-pointing mode, which corresponds to the usual Earth-pointing satellite attitude, and inertial mode (also denoted by IRF = Inertial reference frame) (see Fig. 3).

Fig. 3
figure 3

The two possible configurations: nadir-pointing mode (a) and inertial mode (b)

Full size image

Data were simulated for two configurations of the gradiometer, namely: a single-arm gradiometer, and a double-arm gradiometer. For the single-arm case, gradients were simulated for each of the three possible directions: Txx or Tyy or Tzz; for the double-arm case, gradients were simulated for Txx and Tzz and for Tyy and Tzz. Note that in order to correctly compare two simulations, it was decided that the number of atoms in the cloud had to be kept unchanged in both cases (single-arm and double-arm gradiometer). Obviously, increasing the overall number of atoms would have generally implied a better performance for the double-arm instrument.

For all the simulated observation scenarios, a Monte Carlo (MC) approach was applied in order to compute the estimation errors of the results of the numerical simulations (grid errors and coefficient errors). The MC strategy consisted in simulating sample sets of spherical harmonic model coefficients \(T_{\ell m}\) and then computing error variances and covariances of these solutions. Twenty samples were generated for each scenario (Migliaccio et al., 2009), using the coefficient variances of the reference model EIGEN_6C4, which means that the sample data were equal to EIGEN_6C4 on average, but differ from it by a variation, related to the EIGEN_6C4 error variances. The error r.m.s. of the estimated coefficients \(\hat{T}_{\ell m}\) was thus computed as:

$$\hat{\sigma }_{\ell m} = \sqrt {\frac{1}{20}\mathop \sum \limits_{i}^{20} \left( {\hat{T}_{\ell m}^{i} - T_{\ell m}^{i} } \right)^{2} } \quad \forall { }\ell ,m.$$
(1)

A higher number of samples would be recommendable to get more accurate error estimates; however, here the choice was driven by the need of limiting the computational burden.

2.4 The Cold Atom Gradiometer Noise and Wiener Deconvolution Filter

As described in Sect. 2.2, the first step of the space-wise approach consists in applying a Wiener filter to the data along the orbit to reduce the time correlated noise of the gradiometer. In the case of the MOCASS cold atom gradiometer, time correlation arises also from the atom interferometer integrator or transfer function (Tino & Kasevich, 2014). This means that the MOCASS gradiometer does not provide point-wise observations and, therefore, the Wiener filter was generalized to a Wiener deconvolution filter.

The basic assumption of this filtering technique is that we have a flow of data consisting of the observation signal and noise, which are both stationary processes (with their own respective covariances) and independent from each other. The observation equation in the time domain can be written as:

$$y_{o} \left( t \right) = h\left( t \right)*y\left( t \right) + \nu \left( t \right)$$
(2)

where \(y_{o} \left( t \right)\) is the vector of the observed values; \(h\left( t \right)\) is the atom interferometer integrator; \(y\left( t \right)\) is the signal and \(\nu \left( t \right)\) is the observation noise. The symbol \(*\) stands for the convolution operator.

In the MOCASS numerical simulations, the cold atom interferometer integrator parameters were taken into account (see Fig. 4 and Table 1). More technical details about the interferometer integrator can be found in the papers by Tino and Kasevich (2014) and Migliaccio et al. (2019).

Fig. 4
figure 4

The MOCASS interferometer integrator shape in the frequency (a) and time domain (b)

Full size image
Table 1 The cold atom interferometer parameters
Full size table

The optimal filtering of the data along the orbit was obtained according to the Wiener-Kolmogorov theory (Papoulis, 1984; Sansò, 1986). In the case of GOCE, a Wiener filter was applied (Reguzzoni, 2003), while in the case of MOCASS, a Wiener deconvolution filter was implemented to account for the cold atom interferometer integrator. The MOCASS filter is therefore given by:

$$W\left( f \right) = \frac{{H\left( f \right)S_{y} \left( f \right)}}{{H\left( f \right)^{2} S_{y} \left( f \right) + S_{\nu } \left( f \right)}}$$
(3)

where \(H\left( f \right)\) is the Fourier transform of the interferometer integrator; \(S_{y} \left( f \right)\) is the signal PSD and \(S_{\nu } \left( f \right)\) is the noise PSD. Figure 5 represents the shape of the resulting Wiener deconvolution filter \(W\left( f \right)\) for different observation signals, considering the nadir-pointing mode and high orbit scenario in the case of a single-arm gradiometer. Once the Fourier transform of the observations \(y_{o} \left( t \right)\) is computed, the Wiener deconvolution filter is applied in a computationally efficient way in the frequency domain.

Fig. 5
figure 5

The Wiener deconvolution filter for the MOCASS cold atom gradiometer

Full size image

In the MOCASS numerical simulations, the shape of the signal and noise PSD depended on the simulation scenarios, so different PSD curves were used for the different cases. Examples of the analytical spectra of the noise for the case of the nadir-pointing mode are shown in Fig. 6 for the high orbit (left) and low orbit cases (right). The corresponding curves for the case of the inertial mode are almost identical. As it can be seen, one of distinctive features of the noise PSD of the cold atom gradiometer is the presence of spikes at certain frequency interval (corresponding to about 10–4 Hz and 10–3 Hz). These spikes are related to a residual rotation rate at certain frequencies due to the fact that the satellite attitude cannot be perfectly controlled, so this effect was taken into account in the simulations (Migliaccio et al., 2019).

Fig. 6
figure 6

Analytical spectra of the noise of the MOCASS cold atom gradiometer in nadir-pointing mode in high (left) and low orbit (right): Txx observations (a, b); Tyy observations (c, d); Tzz observations (e, f). In the lower panel, blue and red lines are superimposed

Full size image

In the case of the single-arm gradiometer configuration (for Txx, Tyy and Tzz signals), the noise PSD basically remains the same in each scenario. In the case of the double-arm gradiometer configuration, there is a slight degradation with respect to the single-arm one (the noise is about 1.5 times higher). This is due to the fact that the overall number of cold atoms does not change in the double-arm gradiometer case, but the available cold atoms are split into two directions.

In Fig. 7 the analytical spectrum of the noise is shown against the signal PSD for the cases of Txx, Tyy and Tzz observations (for a single-arm gradiometer), both for the nadir-pointing and inertial mode in high and low orbit.

Fig. 7
figure 7

Analytical spectra of signal and noise of the MOCASS gradiometer (in the case of a single-arm gradiometer): nadir-pointing mode, high orbit (a); nadir-pointing mode, low orbit (b); inertial mode, high orbit (c); inertial mode, low orbit (d)

Full size image

At epoch t, the filtered observation is then represented by the signal plus the filtering error \(e\left( t \right)\):

$$\hat{y}_{o} \left( t \right) = y\left( t \right) + e\left( t \right).$$
(4)

Figure 8a shows the shape of the filtering error covariance function for Txx observations in the case of a single-arm gradiometer, considering different orbit altitudes and pointing modes. As one can see in Fig. 8b for the case of nadir-pointing mode in high orbit, these covariance functions show a quite short correlation that approaches zero within 100 s.

Fig. 8
figure 8

Error covariance function of filtered Txx observations: all scenarios (a) and one of them for a longer time span (b)

Full size image

The parameters of the filtering error covariance function for the case of a single-arm gradiometer and a double-arm gradiometer are reported in Tables 2 and 3, respectively.

Table 2 Formal error of the Wiener deconvolution filter: statistics for the single-arm gradiometer
Full size table
Table 3 Formal error of the Wiener deconvolution filter: statistics for the double-arm gradiometer
Full size table

2.5 Gridding Procedure

After applying the Wiener deconvolution along orbit, a collocation gridding procedure was applied to local patches of data to obtain interpolated values over a boundary sphere Σ at mean satellite altitude. Moreover, a preliminary operation was carried out with the aim of numerically stabilizing the subsequent gridding solution by reducing the signal correlation in such a way that it is consistent with the size of local patches. This consisted in computing a satellite-to-satellite tracking (SST) solution for the spherical harmonic coefficients of the gravity field. In this way, a “MOCASS-only” solution was obtained without combining data with any external a-priori information.

The approach to compute the simulated SST solution for MOCASS was similar to the one used for the analysis of GOCE data (Pail et al., 2011) and consisted in:

  • estimating the anomalous potential from kinematic orbits by the energy integral approach;

  • estimating the spherical harmonic coefficients by least squares adjustment.

The anomalous potential was synthesized along the orbit from the EIGEN_6C4 reference model, adding white noise with 1.35 m2/s2 standard deviation (see Migliaccio et al., 2010). Moreover, since the velocity estimation from kinematic orbit coordinates was based on a moving window of 31-s lag, an under-sampling of 1:31 was applied to comply with the hypothesis that the observation noise is uncorrelated in time. Spherical harmonic coefficients were estimated up to degree 100, also applying a standard Kaula regularization. This maximum spherical harmonic degree was chosen for computational reasons and because the error degree variances of the coefficients estimated from the SST solutions approach the signal degree variances close to degree 100, as shown in Fig. 9.

Fig. 9
figure 9

Error degree variances of the SST solution

Full size image

Coming to the gridding, the local collocation applied for the estimate of interpolated values over the boundary sphere Σ at equiangular grid nodes is represented by:

$$\hat{z}\left( {\vartheta , \lambda } \right) = C_{zy} \left( {C_{yy} + C_{ee} } \right)^{ - 1} \hat{y}_{o}$$
(5)

where \(\hat{y}_{o}\) is the vector of observations forming the local patch of data around the grid node \(P\left( {\vartheta , \lambda } \right)\); \(\hat{z}\left( {\vartheta , \lambda } \right)\) is the predicted functional of the potential at the grid node P; \(C_{zy}\) is the cross-covariance vector between observed and predicted signal; \(C_{yy}\) is the covariance matrix of the signal; \(C_{ee}\) is the covariance matrix of the filtering error; correlations between signal and filtering error are neglected.

In all scenarios, a collocation procedure was applied to the simulated Txx, Tyy and Tzz observations in order to predict values of other functionals of the anomalous potential such as Trr, Tλλ and T on the grid. Basically, these changes of functionals are possible by exploiting the cross-covariance information between observed and predicted quantities. Note that the gridding procedure was implemented by least-squares collocation on local data patches because the full dataset of observations could not be processed as a whole in a unique collocation gridding procedure, due to the huge amount of data.

A particular issue in the gridding procedure regarded the data selection criteria, which is fundamental to ensure the good quality of the results. Some preliminary tests were performed during the study, considering either fixed size equiangular “patches” of data or adaptable “clouds” of data. In the latter case (which was the one applied in the MOCASS study), the idea is to consider clouds of observation points to compute the values of the grid nodes: the closer to the nodes to be estimated, the more observations are used and vice versa, the far from the nodes to be estimated, the less observations are used. In this way the local information is enhanced and possibly all data are used at least once with suitable overlapping clouds. To this aim, the points inside each cloud were randomly chosen by also applying the following conditions: use the largest possible set of input data, and not use the same input data to estimate the values for different grid nodes. The idea of the “clouds” of data is represented in Fig. 10. An equiangular grid with a resolution of 0.5° was devised in the case of MOCASS data processing. Each cloud of data points around the grid nodes is basically divided into three areas with different radius size and different data density.

Fig. 10
figure 10

Example of a cloud of points used as input for the gridding procedure (angular distance measured in degrees); the three areas with different spatial densities are bounded by dashed green lines; red crosses represent grid nodes to be estimated by using this cloud

Full size image

The three areas are:

  • a “full density” (FD) area with the nearest points to the grid nodes to be estimated; no under-sampling is here applied;

  • a “high density” (HD) area with points at intermediate distance to the grid nodes to be estimated;

  • a “low density” (LD) area with the farthest points to the grid nodes to be estimated.

The size of each area depends on the available data, which means that it is not constant for all clouds, but ranges between a minimum and a maximum value, thus adapting to the local density of the data, which in turn mainly depends on the latitude. In the case of MOCASS, Table 4 shows the minimum and maximum size of these areas and the maximum number of points that they could contain. Overall, 8325 clouds of data (containing a maximum of 5000 points) without under-sampling were used (Reguzzoni et al., 2014).

Table 4 Dimensions and number of points in the three areas constituting the “cloud” of points that were used for the data gridding procedure in the MOCASS simulations
Full size table

The global signal covariance function, i.e. the signal degree variances, was adapted to the local signal variances of the “clouds” by applying suitable scale factors. They were computed as the ratio between the local and global variance of Trr signal, even if the observation functional was different. The scale factors for the signal covariance adaptation are shown in Fig. 11; obviously, high values of the scale factors are related to areas where the Trr local signal is much stronger than the global one and, therefore, the collocation filtering is less pronounced.

Fig. 11
figure 11

Scale factors for the local adaptation of the signal covariance function

Full size image

The error r.m.s. of the predicted grid of Trr values is reported in Tables 5 and 6 for the cases of a single-arm and of a double-arm gradiometer.

Table 5 Error r.m.s. of predicted grid of T, Tλλ and Trr values in the case of a single-arm gradiometer; values for optimal case scenarios are highlighted in italics
Full size table
Table 6 Error r.m.s of predicted grid of T, Tλλ and Trr values in the case of a double-arm gradiometer; values for optimal case scenarios are highlighted in italics
Full size table

Tables 5 and 6 show that:

  • for the single-arm gradiometer configuration, observations of Tzz in nadir-pointing mode and observations of Tyy in inertial mode represent the optimal cases, i.e. in these cases the error r.m.s. remains below 0.50 mE;

  • for the double-arm gradiometer configuration, observations of the pair Tyy and Tzz represent the optimal case in both pointing modes, i.e. in this case the error r.m.s. remains below 0.55 mE; we recall that in the case of the double-arm gradiometer, the available cold atoms are split into two directions, and this is why the overall accuracy is not improved with respect to the single-arm gradiometer.

For the sake of comparison, the error r.m.s. of the predicted grids of Trr and Tλλ, computed considering Tzz observations in nadir-pointing mode and high orbit, are displayed in Figs. 12 and 13, respectively. As it can be seen, the Tλλ component is less sensitive to the polar gaps effect than Trr. Note that Tλλ is the second derivative of the anomalous potential T with respect to the longitude. Since the longitude is an angle, the units are not changed, and the result is not formally a “gravity gradient”.

Fig. 12
figure 12

Trr grid error r.m.s., including the polar cap areas (a) and excluding the polar cap areas in order to enhance the visibility of the error (b) (note the color scale in the maps)

Full size image
Fig. 13
figure 13

Tλλ grid error r.m.s., including the polar cap areas (a) and excluding the polar cap areas (b) (note the color scale in the maps)

Full size image

Finally, Fig. 14 shows the grid errors along parallels, in terms of r.m.s. of the Trr predicted data, highlighting the difference of the results when predictions in the polar gaps are also considered. In this figure, the simulated observations are Tzz in nadir-pointing mode and high orbit.

Fig. 14
figure 14

Grid r.m.s. of the Trr predicted data along parallels for all latitudes (a) and excluding polar cap areas (b) (note the different scale of the ordinate [mE] in the two plots)

Full size image

3 Spherical Harmonic Analysis

Once the T, Tλλ and Trr grids have been predicted by the local collocation procedure applied to the simulated and filtered data, spherical harmonic coefficients can be computed by numerical integration (Migliaccio et al., 2004a) for each grid, and then combined according to the MC solution errors to obtain a final estimate of the global geopotential model. Figure 15 shows the complete scheme of the MOCASS data simulations and space-wise analysis.

Fig. 15
figure 15

The complete scheme of the MOCASS data simulations and analysis

Full size image

The spherical harmonic analysis is applied to the gridded observations:

$$\hat{z}\left( {\vartheta , \lambda } \right) = z\left( {\vartheta , \lambda } \right) + \eta \left( {\vartheta , \lambda } \right) = \mathop \sum \limits_{\ell } \mathop \sum \limits_{m} T_{\ell m} Y_{\ell m} \left( {\vartheta , \lambda } \right) + \eta \left( {\vartheta , \lambda } \right)$$
(6)

where \(\eta \left( {\vartheta , \lambda } \right)\) is the gridding error. The estimate of the geopotential coefficients is computed by the following equation:

$$\hat{T}_{\ell m} = \user2{ }\frac{1}{{4\pi a_{\ell } }}\mathop \int \limits_{\Sigma }^{{}} \hat{z}\left( {\vartheta , \lambda } \right)Y_{\ell m} \left( {\vartheta ,\lambda } \right)d\sigma \approx \frac{1}{{4\pi a_{\ell } }}\mathop \sum \limits_{i} \mathop \sum \limits_{j} \hat{z}\left( {\vartheta_{i} , \lambda_{j} } \right)Y_{\ell m} \left( {\vartheta_{i} , \lambda_{j} } \right)\Delta \sigma_{ij}$$
(7)

where \(T_{\ell m}\) is the spherical harmonic coefficient of degree \(\ell\) and order m and \(Y_{\ell m} \left( {\vartheta ,\lambda } \right)\) is the corresponding spherical harmonic function. The parameter \(a_{\ell }\) depends on the type of data. For example, in the case of the second radial derivative of the potential:

$$a_{\ell } = \frac{GM}{{R^{3} }}\left( {\ell + 1} \right)\left( {\ell + 2} \right)\left( \frac{R}{r} \right)^{\ell + 3}$$
(8)

where the constant GM is the product between the gravitational constant G and the mass of the Earth M; R is the reference radius for the Earth and r is the radius of the spherical grid, which in our case is the mean orbital radius Rsat.

4 Gravity Field Recovery Results

Considering 2-month solutions and excluding near-zonal coefficients that are strongly degraded by polar gaps (Sneeuw & van Gelderen, 1997), the results of the MOCASS numerical simulations are shown in Figs. 16 and 17 in terms of error degree variances, for the high orbit case and for the low orbit case, respectively. In these figures, the MOCASS results are presented against those of GOCE and GRACE based on the same amount of data, i.e. considering the error of the GOCE_TIM_R1 model (Pail et al., 2011) and the error of the average of two GRACE monthly solutions (Mayer-Gürr et al., 2016), corresponding to the same data period of the GOCE_TIM_R1 model. The improvement of MOCASS is clearly visible. In fact, in the “optimal” scenarios MOCASS observations allow for an improvement with respect to GOCE at all spherical harmonic degrees and with respect to GRACE at higher harmonic degrees (approximately above 35–40). At lower harmonic degrees, the GRACE estimates remain better than the MOCASS ones.

Fig. 16
figure 16

Error degree variances when estimating a global gravity model for a 2-month mission in high orbit: nadir-pointing mode (a) and inertial mode (b)

Full size image
Fig. 17
figure 17

Error degree variances when estimating a global gravity model for a 2-month mission in low orbit: nadir-pointing mode (a) and inertial mode (b)

Full size image

It has to be underlined that this comparison is a general indication of how well the MOCASS mission would perform with respect to GOCE and GRACE; in fact, the quality of real solutions such as GOCE_TIM_R1 and ITSG-Grace2016 also depends on other factors, e.g., the processing method and the version of the model, which have not been considered in the present study.

For further comparison, the r.m.s. of individual spherical harmonic coefficients estimated in the optimal scenarios from data of a single-arm gradiometer are displayed in logarithmic scale in Fig. 18 for the case of high orbit. The spherical harmonic order is displayed along the horizontal axis, whereas the degree is along the vertical axis. A negative value of the order indicates the sine coefficients, a positive value the cosine ones.

Fig. 18
figure 18

The r.m.s. of the estimated spherical harmonic coefficients in optimal scenarios: Tzz observation in nadir-pointing mode (a) and Tyy observation in inertial mode (b)

Full size image

The results can also be considered in terms of gravity anomaly (Δg) cumulative error at ground level, again for the 2-month solutions. Figure 19 shows the performances of MOCASS for the low orbit case (which of course delivers more informative observations on the gravity field) both in nadir-pointing mode and in inertial mode. The “optimal” results are obtained from observations of the Tzz component in nadir-pointing mode, with an error of 3 mGal at harmonic degree 250, and from observations of the Tyy component in inertial mode, with an error of about 3.5 mGal at harmonic degree 250. The cumulative error was also computed for Trr at satellite altitude and is shown in Fig. 20. Again, the comparison was extended to GOCE and GRACE results, considering the same models used in Figs. 16 and 17. Of course, if the aim is to obtain a static gravity field model, the quality of the solution improves by exploiting a much longer observation period, at least of the duration of 1 year.

Fig. 19
figure 19

Cumulative gravity anomaly error Δg at ground level for a 2-month mission in low orbit: nadir-pointing mode (a) and inertial mode (b)

Full size image
Fig. 20
figure 20

Cumulative error for Trr at satellite altitude for a 2-month mission in low orbit: nadir-pointing mode (a) and inertial mode (b)

Full size image

Figure 21 shows the expected improvement in the estimation of Δg at ground level and Trr at satellite altitude for different mission periods (1-, 2- and 5-year data), considering Tzz observations from a nadir-pointing mission in high orbit. These curves have been obtained by properly rescaling the corresponding ones in Figs. 19 and 20, assuming independent 2-month solutions. For the sake of comparison, the curves for the GOCE and GRACE missions represent results based on 5-year data, obtained from GOCE_TIM_R5 (Brockmann et al., 2014) and ITSG-Grace2014k (Mayer-Gürr et al., 2014) models, respectively.

Fig. 21
figure 21

Cumulative error for MOCASS 1-, 2- and 5-year mission, in terms of Δg at ground level (a) and Trr at satellite altitude (b), considering Tzz observations from a nadir-pointing mission in high orbit

Full size image

Further tests were also performed to assess the level of accuracy attainable with MOCASS if the aim is to estimate the time-variable gravity signal; the results are represented in Fig. 22, which shows the cumulative errors of gravity anomalies Δg at ground level and Trr at satellite altitude. The curves are obtained by propagating the MOCASS monthly solution error to the estimation of a possible linear trend, for a 1-, 2- and 5-year mission, again considering Tzz observations from a nadir-pointing mission in high orbit.

Fig. 22
figure 22

Cumulative errors of gravity anomalies Δg (a) and of Trr (b); linear trend of 1-month solutions, considering Tzz observations from a nadir-pointing mission in high orbit

Full size image

In Fig. 22, the GOCE and GRACE curves represent the linear trend of monthly solutions for a 5-year mission, rescaling the whole mission error of the GOCE_TIM_R5 and ITSG-Grace2014k models on 1 month only. As it can be seen from Fig. 22a, in this case the MOCASS commission error for gravity anomalies is below 1 mGal/month at harmonic degree 300 (and below 0.6 mGal/month at degree 250) even for a 1-year mission; for a 5-year mission it remains below 0.1 mGal/month at harmonic degree 300 (and below 0.05 mGal/month at degree 250).

Figure 23a shows another comparison between MOCASS and the GRACE and GOCE missions based on 2-month data in terms of the error degree variances, including and excluding the effect of polar gaps, i.e. including or excluding error variances of near-zonal coefficients in the computation. As one can see, if polar gaps are included (+ PG), MOCASS gives worse estimates than GRACE, but slightly better ones than GOCE at higher harmonic degrees; however, if polar gaps are excluded (− PG), MOCASS results are better than GOCE at all harmonic degrees, but worse than GRACE at lower harmonic degrees. Therefore, MOCASS estimates at lower harmonic degrees are not satisfactory to detect the time-variable gravity field for geophysical applications; in fact, as an example, in Fig. 23a the MOCASS error is seen to be above the time-variable gravity induced by continental hydrology processes (Iran Pour et al., 2015). Furthermore, a comparison was made with the GOCO05S combined model (Mayer-Gürr et al., 2015), obtained from both GOCE and GRACE data. MOCASS estimates are better than GOCO05S at high harmonic degrees (e.g., over degree 50 for a 5-year mission), as displayed in Fig. 23b.

Fig. 23
figure 23

Error degree variances when estimating a global gravity model: comparison of MOCASS with GRACE and GOCE (a) and comparison of MOCASS with the combined GOCO05S model (b)

Full size image

Finally, Table 7 summarizes the results for the “optimal” MOCASS solutions (in terms of degree variances) for the static and time-variable gravity field, for a mission in high orbit; the GOCE TIM and GRACE solutions corresponding to the MOCASS static and time-variable gravity field solutions are also shown in the table, for the sake of comparison. In the case of static gravity field, the maximum harmonic degree \(\ell\) = 200 is considered, while in the case of time-variable gravity field the maximum degree \(\ell\) = 45 is taken into account, since the signal power is much lower and a higher accuracy is required for the signal detection, at the cost of a lower spatial resolution.

Table 7 Static and time-variable gravity field estimation: comparison of MOCASS “optimal” solutions with GOCE TIM and GRACE solution in terms of Δg at ground level; note that GOCE (TIM_R5) and GRACE (ITSG-Grace2014k) solutions have been rescaled for the corresponding number of months
Full size table

5 Discussion and Conclusions

In the frame of the MOCASS mission proposal, a research was conducted to investigate the impact of a cold atom gradiometer on board a low Earth orbiter acquiring data for the estimation of gravity field models. Several mission scenarios were numerically simulated and the corresponding gravity field models were estimated in order to test the mission performances, using orbit data of about 2 months of the GOCE mission as a “reference orbit”. Important aspects that were considered are:

  • the type of data which could be measured by the MOCASS payload (e.g., second spatial derivatives of the potential Txx, Tyy, Tzz);

  • the characteristics of the MOCASS payload (the cold atom gradiometer) and its noise PSD, from which the along-orbit filtering procedure was defined;

  • the local gridding procedure required to manage the large amount of data by collocation, in turn requiring a smart data partitioning as well as a local estimation of the signal covariance.

The simulations were set up under different configurations of the satellite orbit and pointing mode of the instrument: namely high (259 km) or low (239 km) altitude, and nadir-pointing or inertial mode. The simulated data were Txx, Tyy and Tzz observations for the case of a single-arm gradiometer and a pair of observations (Txx and Tzz) and (Tyy and Tzz) for the case of a double-arm gradiometer.

In terms of reconstruction of a global gravity model, considering all the simulations, the results showed that:

  • in the case of the single-arm gradiometer configuration, Tzz observations in nadir-pointing mode and Tyy observations in inertial mode provided the best accuracy level;

  • in the case of the double-arm gradiometer configuration, Tyy and Tzz observations (in both operation modes) provided a good accuracy level;

  • the results of the simulations of Tzz observations in nadir-pointing mode were also directly compared with GOCE and GRACE results, showing that the MOCASS mission profile could significantly improve the results of the GOCE mission over all the harmonic spectrum, especially at high degrees, with a commission error of about of 1.4 mGal at degree 300 (0.9 mGal at degree 250) for a 5-year mission, but would still be weaker than GRACE at low harmonic degrees, thus limiting its applicability for the time-variable gravity field investigations.