Assessment of gravity field recovery from a quantum satellite mission with atomic clocks and cold atom gradiometers

As a starting point, the observation equation of the potential differences, namely the observable coming from the optical clocks, can be written as:

$$\begin{equation}\Delta {T_0}\left( t \right) = T\left( {\vartheta _t^{\left( i \right)},\lambda _t^{\left( i \right)},r_t^{\left( i \right)}} \right) - T\left( {\vartheta _t^{\left( j \right)},\lambda _t^{\left( j \right)},r_t^{\left( j \right)}} \right) + \varepsilon \left( t \right)\end{equation} \tag{ 5 }$$

where $T$ is the anomalous potential, $\vartheta _t^{\left( \cdot \right)},\lambda _t^{\left( \cdot \right)},r_t^{\left( \cdot \right)}$ are the spherical coordinates of satellite $i$ or $j$ along their orbits at a specific time $t$, and $\varepsilon$ is the atomic clock noise. Given the observation equation (5), we recall that the anomalous potential $T$ can be written in terms of spherical harmonic coefficients as

$$\begin{equation}T\left( {\vartheta _t^{\left( \cdot \right)},\lambda _t^{\left( \cdot \right)},r_t^{\left( \cdot \right)}} \right) = \frac{{GM}}{R}\mathop \sum \limits_\ell {\left( {\frac{R}{{r_t^{\left( \cdot \right)}}}} \right)^{\ell + 1}}\mathop \sum \limits_m {T_{\ell m}}{Y_{\ell m}}\left( {\vartheta _t^{\left( \cdot \right)},\lambda _t^{\left( \cdot \right)}} \right)\end{equation} \tag{ 6 }$$

where $G$ is the gravitational constant, $M$ is the mass of the Earth, $R$ is the reference radius of the Earth, ${Y_{\ell m}}$ is the spherical harmonic function of degree $\ell$ and order $m$, and ${T_{\ell m}}$ is the harmonic coefficient. The overall set of coefficients in equation (6) can be estimated by a least squares adjustment. To this aim, the deterministic model of the system can be written as

$$\begin{equation}{\boldsymbol{\Delta }}{\boldsymbol{{T}}} = {\mathbf{A}}\,{\boldsymbol{{T}}}\end{equation} \tag{ 7 }$$

where ${\boldsymbol{\Delta }}{\boldsymbol{{T}}}$ is the vector of the observables, namely the potential differences, ${\boldsymbol{{T}}}$ the vector of the unknown ${T_{\ell m}}$ spherical harmonic coefficients, and ${\mathbf{A}}$ the design matrix where the number of rows equals the number of observations and the number of columns equals the number of unknowns, depending on the maximum spherical harmonic degree $L$. The single element of matrix ${\mathbf{A}}$ related to the observation at time $t$ for the coefficient of degree $\ell$ and order $m$, considering the two satellites $i$ and $j$ can be written as

$$\begin{equation}\begin{array}{*{20}{r}} {{a_{t,\ell m}}}&{ = \dfrac{{GM}}{R}\left[ {{{\left( {\dfrac{R}{{r_t^{\left( i \right)}}}} \right)}^{\ell + 1}}{Y_{\ell m}}\left( {\vartheta _t^{\left( i \right)},\lambda _t^{\left( i \right)}} \right) - {{\left( {\dfrac{R}{{r_t^{\left( j \right)}}}} \right)}^{\ell + 1}}{Y_{\ell m}}\left( {\vartheta _t^{\left( j \right)},\lambda _t^{\left( j \right)}} \right)} \right]} \end{array}.\end{equation} \tag{ 8 }$$

Regarding the stochastic model, we considered white observation noise, with zero mean and known variance $\sigma _\varepsilon ^2$:

$$\begin{equation}E\left[ {\varepsilon \left( t \right)} \right] = 0\ {\text{and}}\ {{\mathbf{C}}_{\varepsilon \varepsilon }} = \sigma _\varepsilon ^2{\mathbf{I}}.\end{equation} \tag{ 9 }$$

The variance is defined according to the assumption presented in section 2. Since potential difference instead of potential is observed, the problem is ill-posed. This requires to apply a proper regularization procedure with the additional advantage of numerically stabilizing the estimates of high degree coefficients. For this purpose, either the well-known Kaula rule (Kaula 1966) or an a-prior degree variance model can be applied, introducing a regularization matrix ${\mathbf{R}}$ equal to

$$\begin{equation}{\mathbf{R}} = \left[ {\begin{array}{*{20}{c}} \ddots &0&0&0 \\ 0&{\sigma _\ell ^2}&0&0 \\ 0&0& \ddots &0 \\ 0&0&0& \ddots \end{array}} \right]\end{equation} \tag{ 10 }$$

where $\sigma _\ell ^2$ is the degree variance of the spherical harmonic coefficient ${T_{\ell m}}$.

The least-squares solution, directly estimating the spherical harmonic coefficients, is given by

$$\begin{equation}\hat T = \left[ \begin{gathered} \;\; \vdots \\ {{\hat T}_{\ell m}} \\ \;\; \vdots \\ \end{gathered} \right] = \left( {{\mathbf{A}^{\text{T}}}\mathbf{C}_{\varepsilon \varepsilon }^{-1}\mathbf{A} {\text{ + }}{\mathbf{R}^{ - 1}}} \right){\mathbf{A}^{\text{T}}}\mathbf{C}_{\varepsilon \varepsilon }^{-1}{\boldsymbol{\Delta} {\boldsymbol{T}_0}}\end{equation} \tag{ 11 }$$

where ${\mathbf{R}}$ is the regularization matrix of equation (10) and ${\boldsymbol{\Delta }}{{\boldsymbol{{T}}}_0}$ is the vector containing all the available observations. Then, the error covariance matrix of the estimated coefficients can be expressed as

$$\begin{equation}{{\mathbf{C}}_{\hat T\hat T}} = {\left( {{{\mathbf{A}}^{\text{T}}}{\mathbf{C}}_{\varepsilon \varepsilon }^{ - 1}{\mathbf{A}} + {{\mathbf{R}}^{ - 1}}} \right)^{ - 1}}.\end{equation} \tag{ 12 }$$

Wiener filter or deconvolution (Papoulis 1981, Reguzzoni 2003, Albertella et al 2004) is a necessary step in the space-wise data processing since it is applied to reduce the time-correlated noise of the observations before the gridding. A short error correlation in time (and therefore in space) is required for gridding data on local areas, while a low error variance leads to more accurate collocation estimates. The use of either a Wiener filter or a Wiener deconvolution depends on the fact that the observations are either point-wise (as for the electrostatic gradiometers) or the outcome of an integration in time (as it happens for the quantum gradiometers based on CAI).

Based on the gradient observations proposed for MOCAST+, it is necessary to devise a deconvolution scheme accordingly. In particular, in the time domain the observation equation of the second derivatives of the gravity field is written as

$$\begin{equation}{y_o}\left( t \right) = h\left( t \right)*y\left( t \right) + \nu \left( t \right)\end{equation} \tag{ 13 }$$

where ${y_o}\left( t \right)$ is the vector containing the observed second derivatives (given the direction of the arm of the gradiometer), $h\left( t \right)$ is the transfer function of the instrument (Migliaccio et al 2019), $y\left( t \right)$ is the noiseless gravity gradient signal, the $*$ operator stands for the time convolution, and $\nu \left( t \right)$ is the gradiometer measurement noise.

Starting from the analytical signal and noise ASDs presented in section 2, and also introducing the cold atom gradiometer transfer function (Migliaccio et al 2019), the Wiener deconvolution operator $W\left( f \right)$ can be derived by minimizing the mean square estimation error, thus obtaining

$$\begin{equation}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)}}\end{equation} \tag{ 14 }$$

where $H\left( f \right)$ is the Fourier transform of the transfer function, ${S_y}\left( f \right)$ is the power spectral density (PSD) of the signal $y\left( t \right)$ and ${S_\nu }\left( f \right)$ is the PSD of the noise $\nu \left( t \right)$. Note that the PSD is computed as the square of the ASD. The filtered signal $\hat y\left( t \right)$ can be retrieved by applying the deconvolution operator $W\left( f \right)$ to the Fourier transform of the observed second derivatives ${y_o}\left( t \right)$, together with its estimation error. The PSD of this error can be obtained as

$$\begin{equation}{S_{\hat e}}\left( f \right) = \frac{{{S_y}\left( f \right){S_\nu }\left( f \right)}}{{H{{\left( f \right)}^2}{S_y}\left( f \right) + {S_\nu }\left( f \right)}}\end{equation} \tag{ 15 }$$

from which the corresponding covariance function can be computed by the inverse Fourier transform. Given this covariance function, the ${{\mathbf{C}}_{\hat e\hat e}}$ covariance matrix of the estimation error of the filtered signal $\hat y\left( t \right)$ can be computed and used in the collocation gridding.

The operator of equation (14) is designed by introducing the hypothesis of having a stationary process in time, with independent observation noise for instruments mounted on different satellites. Indeed, the signal observed from the different satellites is correlated, since they all fly under the action of the same gravity field. In principle this correlation may be exploited to jointly filter/deconvolve observations from all gradiometers; however, it was decided to introduce this information in the gridding procedure only, because this correlation is spatially dependent and therefore the signal covariance and cross-covariance functions are only approximately stationary in time. In other words, the Wiener deconvolution is only a pre-processing step with some approximated hypotheses, while the core of the space-wise approach is the local gridding by collocation, where the spatial signal correlation and the time noise correlation can be both properly modelled.

The Wiener deconvolution is followed by a collocation procedure that allows to obtain gridded data starting from the observed gravity gradients. For the proposed mission, ${T_{xx}}$, ${T_{yy}}$, and ${T_{zz}}$ were considered as the observed gradients, where $x$ is the along-orbit direction, $y$ is the out-of-orbital-plane direction, and $z$ is the radial direction. Note that if the orbit is not exactly circular, the $x$ direction is 'almost' along the orbit in order to maintain the $z$-axis along the radial direction. The functionals to be predicted by the gridding procedure were the potential $T$, the second derivative of the potential along the radial direction ${T_{rr}}$, and the second derivative of the potential with respect to the longitude ${T_{\lambda \lambda }}$. Note that the last quantity is not formally a gravity gradient because it is a second derivative with respect to an angle and has the same measurement unit of the potential $T$. As a starting point, the observation equation corresponding to the position at epoch $\hat y\left( {{\vartheta _t},{\lambda _t},{r_t}} \right) = y\left( {{\vartheta _t},{\lambda _t},{r_t}} \right) + {\hat e_t}$ is written as

$$\begin{equation}\hat y\left( {{\vartheta _t},{\lambda _t},{r_t}} \right) = y\left( {{\vartheta _t},{\lambda _t},{r_t}} \right) + {\hat e_t}\end{equation} \tag{ 16 }$$

where $\left( {{\vartheta _t},{\lambda _t},{r_t}} \right)$ are the coordinates at time $t$, the vector $\hat y$ contains the Wiener filtered gradients, and ${\hat e_t}$ is the estimation error of the Wiener deconvolution, see equations (14) and (15).

The gridding procedure is done by collocation over local patches of data, since a global collocation procedure would be numerically unfeasible (Migliaccio et al 2006, Reguzzoni et al 2014). Moreover, with local collocation it is possible to avoid numerical drawbacks due either to considering a too large amount of data at once, or to including observation points which are too close to each other. The estimated grid is defined on a spherical boundary at a constant altitude equal to the mean satellite altitude (or at least included between the lower and higher orbit altitudes, in case of many satellites), and each grid node is computed by considering the neighbouring data. Therefore, one of the most crucial parts of the gridding procedure is the determination of the local patches of data based on the position of the nodes. To enhance the local information, more observations are used in the proximity of the grid node to be estimated, also applying the condition that the same set of points must not be used to estimate different grid nodes.

The partitioning of the set of data points to be considered for the estimation of grid nodes is done by dividing the point patch into three areas with different size and point density. The size of the patch depends on the correlation length of the covariance function of the data, while the size of each area varies depending on the amount of the available data points, implying a possible data under-sampling. A representation of this 'adaptable' patch of data points can be seen in figure 4.

Zoom In Zoom Out Reset image size

Before proceeding with the collocation step, the long-wavelength effects must be removed from equation (16), so the clock-only solution computed by equation (11) is subtracted from the original observation equation, obtaining

$$\begin{equation}\delta \hat y\left( {{\vartheta _t},{\lambda _t},{r_t}} \right) = \hat y\left( {{\vartheta _t},{\lambda _t},{r_t}} \right) - \overline y \left( {{\vartheta _t},{\lambda _t},{r_t}} \right)\end{equation} \tag{ 17 }$$

where $\overline y \left( {{\vartheta _t},{\lambda _t},{r_t}} \right)$ represents the gravity signal computed by the synthesis of the ${\hat T_{\ell m}}$ spherical harmonic coefficients estimated by the least squares solution of equation (11). Recalling both equation (16) and the following formula

$$\begin{equation}y\left( {{\vartheta _t},{\lambda _t},{r_t}} \right) = \overline y \left( {{\vartheta _t},{\lambda _t},{r_t}} \right) + \delta y\left( {{\vartheta _t},{\lambda _t},{r_t}} \right)\end{equation} \tag{ 18 }$$

where $\delta y\left( {{\vartheta _t},{\lambda _t},{r_t}} \right)$ is the error of the gradients synthesised from the clock-only solution, namely the missing signal to be retrieved by the gridding procedure, the observation equation can be written as

$$\begin{equation}\delta \hat y\left( {{\vartheta _t},{\lambda _t},{r_t}} \right) = \delta y\left( {{\vartheta _t},{\lambda _t},{r_t}} \right) + {\hat e_t}.\end{equation} \tag{ 19 }$$

From here, the collocation estimate can be derived as

$$\begin{equation}\hat z\left( {\vartheta ,\lambda } \right) = {{\mathbf{C}}_{z\delta y}}{\left( {{{\mathbf{C}}_{\delta y\delta y}} + {{\mathbf{C}}_{\hat e\hat e}}} \right)^{ - 1}}{\boldsymbol{\delta}} {\boldsymbol{\hat y}}_{o}\end{equation} \tag{ 20 }$$

where $\hat z\left( {\vartheta ,\lambda } \right)$ is the functional of the gravitational field (i.e. $T$, ${T_{rr}}$, or ${T_{\lambda \lambda }}$) to be estimated at the grid node with coordinates $\left( {\vartheta ,\lambda } \right)$ starting from the given set of observations ${\boldsymbol{{\delta }}}{{\boldsymbol{{\hat y}}}_o}$, ${{\mathbf{C}}_{\delta y\delta y}}$ is the covariance matrix of the observed functionals in the neighbourhood of the grid node, ${{\mathbf{C}}_{z\delta y}}$ is the cross-covariance matrix between the observed functional and the one to be estimated, and ${{\mathbf{C}}_{\hat e\hat e}}$ is the covariance matrix of the Wiener deconvolution estimation error. Note that ${{\mathbf{C}}_{\delta y\delta y}}$ and ${{\mathbf{C}}_{z\delta y}}$ can be computed from equation (12) and locally adapted based on the data inside each point patch, while ${{\mathbf{C}}_{\hat e\hat e}}$ is computed according to equation (15). Also note that, whereas the signal $y$ (and therefore $\delta y$) is uncorrelated from the gradiometer noise $\nu \left( t \right)$, the Wiener deconvolution introduces a correlation between the signal and the estimation error $\hat e$. This correlation can be properly considered by iterating the general scheme of the space-wise approach (Reguzzoni and Tselfes 2009), but here it has been neglected to speed up the computations, given the need of performing a large number of simulations. This approximation is expected to mainly affect the estimation of the highest spherical harmonic degrees.

Considering the observation equations (16) and (19) and rewriting them for the estimated grid values, we obtain

$$\begin{equation}\hat z\left( {\vartheta ,\lambda } \right) = z\left( {\vartheta ,\lambda } \right) + \eta \left( {\vartheta ,\lambda } \right) = \mathop \sum \limits_\ell \mathop \sum \limits_m {a_{\ell m}}\,{T_{\ell m}}{Y_{\ell m}}\left( {\vartheta ,\lambda } \right) + \eta \left( {\vartheta ,\lambda } \right)\end{equation} \tag{ 21 }$$

where $z\left( {\vartheta ,\lambda } \right)$ is the estimated functional (i.e. $T$, ${T_{rr}}$, or ${T_{\lambda \lambda }}$) and $\eta \left( {\vartheta ,\lambda } \right)$ is the gridding error. The gridding procedure to compute each of the estimated functional directly includes the combination of the observations from all satellites, meaning that all the observed gradients falling inside one patch are directly combined through equation (20) to estimate grids of $T$, ${T_{rr}}$, or ${T_{\lambda \lambda }}$ values.

For each of the grids, the spherical harmonic coefficients are derived by numerical integration, computed over the spherical boundary, based on quadrature equations (Colombo 1981)

$$\begin{equation}\hat T_{\ell m}^{\,z} = \frac{1}{{4\pi a_{\ell m}^z }}\int_\Sigma{\hat z\left( {\vartheta ,\lambda } \right){Y_{\ell m}}\left( {\vartheta ,\lambda } \right){\text{d}}\sigma \approx \frac{1}{{4\pi a_{\ell m}^z}}\sum\limits_i {\sum\limits_j {\hat z} } \left( {{\vartheta _i},{\lambda _j}} \right){Y_{\ell m}}\left( {{\vartheta _i},{\lambda _j}} \right)\Delta {\sigma _{ij}}} \end{equation} \tag{ 22 }$$

where the indexes $i$ and $j$ are used to identify the grid nodes, ${\Delta }{\sigma _{ij}}$ is the size of the equiangular grid cell around each the grid node $\left( {i,j} \right)$, and $a_{\ell m}^z$ is the (transfer) function that depends on the type of $z$ functional, which is equal to

$$\begin{equation}a_{\ell m}^T = \frac{{GM}}{R}{\left( {\frac{R}{r}} \right)^{\left( {\ell + 1} \right)}}\end{equation} \tag{ 23 }$$

$$\begin{equation}a_{\ell m}^{{T_{rr}}} = \frac{{GM}}{{{R^3}}}\left( {\ell + 1} \right)\left( {\ell + 2} \right){\left( {\frac{R}{r}} \right)^{\left( {\ell + 3} \right)}}\end{equation} \tag{ 24 }$$

$$\begin{equation}a_{\ell m}^{{T_{\lambda \lambda }}} = - \frac{{GM}}{R}{m^2}{\left( {\frac{R}{r}} \right)^{\left( {\ell + 1} \right)}}.\end{equation} \tag{ 25 }$$

The spherical harmonic coefficients computed from the three estimated grids are then combined to obtain a unique solution. This combination can be performed by weighting the coefficients from different grids according to their estimation errors computed by a Monte Carlo approach (Migliaccio et al 2009). According to this approach, the kth sample of the spherical harmonic coefficients $T_{\ell m}^k$ was drawn from a reference global gravity model (considering its error degree variances) and used to simulate observations along the orbit together with a random sample of the instrumental noise. Then, the simulated observations entered the space-wise computation scheme (figure 3) and grid values ${\hat z^k}$ with the corresponding coefficients $\hat T_{\ell m}^{z,k}$ were estimated for each sample $k$, so that the empirical error rms can be computed as

$$\begin{equation}\hat \sigma _{\ell m}^z = \sqrt {\frac{1}{N}\mathop \sum \limits_{k = 1}^N {{\left( {\hat T_{\ell m}^{\,z,k} - T_{\ell m}^k} \right)}^2}} \,\,\,\forall { }\ell ,m\end{equation} \tag{ 26 }$$

where $N$ represents the overall number of Monte Carlo samples. The coefficient variances computed by equation (26) were used to determine the combination of the spherical harmonic coefficients coming from grids of different functionals as:

$$\begin{equation}\hat T_{\ell m}^k = \frac{{\frac{{\hat T_{\ell m}^{T,k}}}{{{{\left( {\hat \sigma _{\ell m}^T} \right)}^2}}} + \frac{{\hat T_{\ell m}^{{T_{\lambda \lambda }},k}}}{{{{\left( {\hat \sigma _{\ell m}^{{T_{\lambda \lambda }}}} \right)}^2}}} + \frac{{\hat T_{\ell m}^{{T_{rr}},k}}}{{{{\left( {\hat \sigma _{\ell m}^{{T_{rr}}}} \right)}^2}}}}}{{\frac{1}{{{{\left( {\hat \sigma _{\ell m}^T} \right)}^2}}} + \frac{1}{{{{\left( {\hat \sigma _{\ell m}^{{T_{\lambda \lambda }}}} \right)}^2}}} + \frac{1}{{{{\left( {\hat \sigma _{\ell m}^{{T_{rr}}}} \right)}^2}}}}}\forall \ell ,m.\end{equation} \tag{ 27 }$$

Given this combination for each of the $N$ samples, the empirical error rms of the resulting spherical harmonic coefficients can be computed as

$$\begin{equation}\hat \sigma _{\ell m}^{} = \sqrt {\frac{1}{N}\mathop \sum \limits_{k = 1}^N {{\left( {\hat T_{\ell m}^k - T_{\ell m}^k} \right)}^2}} \forall {\text{ }}\ell ,m.\end{equation} \tag{ 28 }$$

Note that the combination in equation (27) is not optimal because covariances and cross-covariances are neglected due to the fact that their computation would require a very high number of Monte Carlo samples.

As for the grids, the error rms at each node can be similarly computed as

$$\begin{equation}{\hat \sigma _z}\left( {{\vartheta _i},{\lambda _j}} \right) = \sqrt {\frac{1}{N}\mathop {{{\mathop \sum \nolimits }}}\limits_{k = 1}^N {{\left[ {{{\hat z}^k}\left( {{\vartheta _i},{\lambda _j}} \right) - {z^k}\left( {{\vartheta _i},{\lambda _j}} \right)} \right]}^2}}\ \forall i,j.\end{equation} \tag{ 29 }$$

For the clock only simulations, three main satellite orbit configurations were considered. They consist of a GRACE-class formation with a pair of satellites on a polar orbit, a Bender-class formation with two pairs of satellites on polar and inclined orbits, and finally a Bender-class formation with three satellites along each orbit.

Each one of these formations was tested for different sampling rates, different inter-satellite distances and, only in the case of the GRACE-class formation, for different satellite altitudes. As for the last point, increasing or lowering the satellite orbit altitude by 120 km with respect to the nominal one shows that the lower orbit scenario, as expected, can provide an improvement at medium-high harmonic degrees (see yellow line in figure 5), while increasing it would have a significant negative effect on the quality of the estimation (see blue line in figure 5). However, lowering the altitude has a limited impact on the low harmonic degrees (the most important ones for time-variable gravity field analysis) and would shorten the mission lifetime. Therefore, this parameter was neglected in the subsequent simulations.

Zoom In Zoom Out Reset image size

Scenarios based on a Bender-class formation allow to improve the results with respect to the GRACE-class formation. Some significant results are reported in figure 6, showing the accuracy of different simulations in terms of error degree variances, computed according to equation (28). Note that the Monte Carlo empirical error variances are basically equivalent to the formal error variances taken from the diagonal of the covariance matrix in equation (12). Considering the Bender-class formation with two satellites per orbit, the quality of the estimated solution can be improved by about two orders of magnitude (in terms of variance) by increasing the inter-satellite distance from 100 km to 1000 km (according to the theoretical limits for this technology, as described in section 2), as it can be seen by comparing the blue and yellow lines in figure 6. Moreover, if the clock observation sampling rate is increased from 0.1 Hz to 1 Hz, this leads to an improvement of about one order of magnitude, as it can be seen by comparing the solid and dashed lines in figure 6. Finally, simulations with an additional third satellite on both the polar and the inclined orbit were performed, in order to increase the spatial resolution and the number of observations. In this case the inter-satellite distance between pairs of satellites was set equal to 1000 km, allowing to slightly improve the results of the previous configurations, as it can be seen in figure 6 by comparing the purple line with the other ones. The three cases presented in figure 6 represent the best-case scenarios for the clock-only solutions, and they are compared with the error degree variances from GRACE. The latter were obtained by the monthly average of ITSG-GRACE 2018 errors (Kvas et al 2019) over a time span of about 14 years. As it can be seen, unfortunately the usage of atomic clocks only is not sufficient to ensure a gravity field solution better than the GRACE one. Only the Bender configuration with six satellites and a clock observation frequency of 1 Hz (dashed purple line in figure 6) could approach the GRACE error level, improving a few very low harmonic degrees.

Zoom In Zoom Out Reset image size

Coming to the main part of the payload, which is the CAI, the simulations were aimed at investigating the impact of the gradiometer arm orientation, which is related to the orbital inclination (polar or inclined) as well as to the different noise PSDs (in turn depending on the level of accuracy which could be achieved with the attitude compensation of the atmospheric drag and other torques, see again figure 1). For this purpose, a period of one month was considered for each simulation; no other gravitational or non-gravitational effects (tidal and non-tidal mass variations) were taken into account.

Before presenting the results of the simulations, it is important to stress the impact of the accuracy of the attitude control sensor on the gravity gradient measurement in the orbital plane (thus involving gradiometers mounted in the along-track $x$-axis and in the radial $z$-axis directions). In fact, angular rotations, through the centrifugal term, put a serious limitation to these measurements. For this reason, compensation of the residual angular rotations around the out-of-orbital-plane $y$-axis direction is needed; the accuracy of this compensation highly influences the gradiometer accuracy on the along-track $x$-axis and radial direction $z$-axis, as already shown in section 2. To understand the impact of changing the gradiometer arm direction ($x$, $y$, $z$), the solutions of a series of scenarios considering observations from a single gradiometer mounted onboard a single satellite (along the polar or inclined orbit) were computed, introducing all possible combinations between gradiometer arm direction and drag compensation level. To allow for a proper comparison, the same long-wavelength model of the gravity field was removed before the gridding procedure for each simulation. For this purpose, the clock-only solution for the Bender-class formation with two pairs of satellites at 100 km inter-satellite distance and 0.1 Hz clock sampling rate (solid blue line in figure 6) was used.

The results for the polar and inclined orbits are shown in figures 7 and 8, respectively. The best solution was obtained by orienting the gradiometer arm along the radial direction ($z$-axis), considering the scenario with a full attitude compensation level. As already discussed in section 2, this is an ideal scenario, therefore the best solution is actually obtained by orienting the gradiometer arm along the out-of-orbital-plane direction ($y$-axis), since in this case the instrumental accuracy is not affected by the accuracy of the attitude control. Note that the inclined orbit is affected by quality degradation due to the data gaps at high latitudes, which lead to a solution degradation between degrees 20 and 120, where the error degree variances are dominated by the higher error of the lower order coefficients. For this reason, to better evaluate the accuracy achievable by exploiting different gradiometer directions, the degree variances were replaced by the squared degree medians of the coefficient error, which are shown in figure 8(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A further comparison of the impact of different gradiometer axis orientations was performed at the level of the single spherical harmonic coefficient, considering only the most realistic case with 1 nrad s⁻¹ attitude control accuracy. The results of this comparison are shown in figure 9. Here, the gradiometer arm direction providing the best estimate, i.e. providing the smallest error variance, is shown for each harmonic degree and order. With this level of attitude control accuracy, the main contribution to the coefficient estimation comes from the gradiometer oriented along the $y$-axis direction, as it can also be observed from figures 7 and 8. This finding is in agreement with previous works on quantum gradiometry, see e.g. Douch et al (2018). Nevertheless, a contribution coming from ${T_{zz}}$ at low orders can be seen on the polar orbit. As for the inclined orbit this contribution is still present, but extremely limited if compared to ${T_{yy}}$.

Zoom In Zoom Out Reset image size

The altitude at which gridded data are predicted is another parameter to be considered in the possible combinations. To check its impact, the following values of the grid altitude were tested: 347 km (mean altitude of the inclined orbit), 371 km (mean altitude of the polar orbit), and 360 km (overall average altitude of the two orbits). However, the results (computed for a single gradiometer on a single satellite) showed a negligible impact of the choice of altitude, therefore in the following simulations the grids were always computed at the average altitude of 360 km.

Given the results of the single-satellite gradiometer simulations, the analysis focused on Bender-class formations with pairs or triplets of in-line satellites, always setting the inter-satellite distance at about 1000 km, thus considering only the best conditions for the clock-only solutions, according to the results presented in figure 6. The comparisons were performed between the best ideal scenario (full attitude compensation) and the more realistic one (1 nrad s⁻¹ attitude compensation). In the last scenario, gradiometers with different atomic species were considered too. Under these conditions, three mission configurations were simulated based on Bender-class formations with two or three satellites per orbit:

• ideal, with full attitude compensation, mounting ${T_{zz}}$ gradiometers on all satellites;
• realistic with Sr gradiometer, with 1 nrad s⁻¹ attitude compensation, mounting ${T_{yy}}$ gradiometers on all satellites apart from a single ${T_{zz}}$ gradiometer on the polar orbit;
• realistic with Rb gradiometer, with 1 nrad s⁻¹ attitude compensation, mounting ${T_{yy}}$ gradiometers on all satellites apart from a single ${T_{zz}}$ gradiometer on the polar orbit.

The gradiometer arm directions were chosen based on the error curves presented in figures 7 and 8, while the choice of mounting a ${T_{zz}}$ gradiometer on a satellite along the polar orbit for realistic scenarios was driven by the results presented in figure 9. A further scenario was devised for the three-satellites-per-orbit configuration, based on one ${T_{yy}}$ gradiometer on board the central satellite of each orbit, in order to minimize the complexity and cost of the mission. It seems reasonable to assume that the central satellite plays the role of master with the gradiometer and two tracking clock systems pointing towards the other two satellites; anyway, the choice of the master does not affect the results of the analysis. The least squares clock-only solution needed to reduce the long-wavelength required by the gridding procedure was computed according to the satellite configuration of each scenario, always considering 1 Hz atomic clock sampling rate.

The results of these simulations are shown in figure 10 (Bender, two satellites per orbit) and figure 11 (Bender, three satellites per orbit), where the empirical error degree variances from Monte Carlo simulations are also compared with the error degree variances from GRACE.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The error curves of figures 10 and 11 show that the Bender-class formation with two or three satellites per orbit allows to estimate spherical harmonic coefficients up to about degree 200 with only one month of observations, irrespectively of the chosen mission configuration. This is quite a remarkable result, when compared with GRACE and GOCE performances. Note that the best solution was obtained when using ideal conditions with a full attitude control (see red line in figures 10 and 11) and represents a sort of theoretical limit for this kind of instruments and orbit configurations. As for the impact of using Sr instead of Rb as atomic species, an improvement by a factor ∼2 in terms of standard deviation was obtained above degree 20, as shown by the comparison of the yellow and purple lines in figures 10 and 11. This improvement is mainly due to the lower level of the Sr gradiometer noise PSD (see figure 1), while the lower degrees are also controlled by the clock-only solution, which is the same for the two scenarios. The configuration with only two ${T_{yy}}$ Sr gradiometers shows an error which, below degree 100, is of the same level of the one with six Rb gradiometers, confirming the possibility of reducing costs and mission complexity without losing too much in terms of accuracy.

A further test was performed to assess the impact of the clocks on the computed solutions. To this aim, the realistic configuration with six satellites (attitude control at 1 nrad s⁻¹ level with five ${T_{yy}}$ and one${ }{T_{zz}}$ gradiometers) was considered introducing a degraded accuracy in the potential difference observations (by a factor 10 in terms of noise variance). Note that the same degradation can be obtained by working at 0.1 Hz with the nominal clock accuracy. The result of this degraded scenario is compared with the one previously computed with the same satellite configuration and the nominal clock accuracy of 0.2 ${{\text{m}}^2}\,{{\text{s}}^{ - 2}}$ (see figure 12). The differences between the two solutions are in the lower degrees, namely below 50, where better measurements from clocks can improve the quality of the overall solution up to one order of magnitude (compare the red and blue solid lines in figure 12). However, the clock-only solution never 'dominates' the error curve at low degrees. Therefore, the improvement could be partially related to the role of the clock-only solution in the space-wise approach, particularly in the data reduction before the gridding affecting the signal covariance modelling (see equations (17)–(20)).

Zoom In Zoom Out Reset image size

For the sake of completeness, the cumulative commission errors for geoid undulations and gravity anomalies at ground level (in spherical approximation) were computed and are shown in figures 13 and 14 for the ideal scenario (${T_{zz}}$ on all the satellites with a full attitude control) and for the realistic scenario (${T_{yy}}$ on all satellites apart from the one with ${T_{zz}}$ considering a 1 nrad s⁻¹ attitude control), either in the case of two or three satellites per orbit in Bender-class formation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, the impact of longer time span solutions was tested for the Bender-class formation in the nominal mission case (one pair of satellites along each orbit at 100 km inter-satellite distance). The error grids of the different time span solutions, computed as the standard deviation of grid errors over 100 Monte Carlo samples, are shown in figure 15 in the case of the second radial derivative ${T_{rr}}$. As it can be seen, increasing the period covered by the solution leads to a more uniform precision, thanks to the denser data distribution. In fact, the one month and two months gravity field solutions still show a visible pattern coming from the orbit ground tracks, while this pattern is strongly reduced in the four-month solution. Additionally, figure 15(d) shows the root mean square error for each point of the grid (for the four-month solution only), thus also showing systematic effects that are not correctly retrieved in the solution and, therefore, are present in all Monte Carlo samples. As expected, some high frequencies are not estimated, mainly depending on the satellite altitude and gradiometer accuracy, as well as on the processing strategy, which was not optimized for the estimation of such high-degree coefficients (i.e. too large clouds of data).

Zoom In Zoom Out Reset image size

The MOCAST+ study was completed by investigating the effects of the non-tidal mass variation on the gravity field generated by the variability in atmosphere (A component), oceans (O component), continental ice-sheets (I component), groundwater (H component), and solid Earth (S component) on the estimation accuracy of monthly solutions, considering different mission scenarios. One of the main problems caused by the time-variable component is the presence of aliasing in the computed solution. The aliasing effect is related to the gravity signal with frequency higher than one month (e.g. from hours to days in atmosphere, oceans, and, considering smaller spatial extent, also in groundwater) and background model errors of these signals. The magnitude of this aliasing effect is usually higher than the noise level of the instruments, which are the CAI for the proposed MOCAST+ mission or the K-band ranging system for GRACE and GRACE-FO missions (Wiese 2011, Wiese et al 2012).

For the purposes of this investigation, three mission configurations were tested, namely:

• a GRACE-class formation, with pair of satellites on a polar orbit, with one ${T_{zz}}$ gradiometer and one ${T_{yy}}$ gradiometer;
• a Bender-class formation with a pair of satellites per orbit, with one ${T_{zz}}$ gradiometer on the polar orbit and three ${T_{yy}}$ gradiometers on the inclined orbit;
• a Bender-class formation with a triplet of satellites per orbit, with one ${T_{zz}}$ gradiometer on the polar orbit and five ${T_{yy}}$ gradiometers on the inclined orbit.

In all cases the clocks were assumed to work at 1.0 Hz sampling rate and the attitude control was assumed to have a 1 nrad s⁻¹ accuracy. For all simulations only non-tidal mass variations were considered by introducing the ESA updated ESM (Dobslaw et al 2015). In particular, the overall effect due to AOHIS was added to the static gravity field of each Monte Carlo sample and reduced with the de-aliasing model (DEAL) provided with the updated ESM. To properly account for model errors in the AO components of DEAL, this term was corrected by the AOerr model (also provided with the updated ESM), thus allowing to simulate a realistic scenario with a consistent dataset (Dobslaw et al 2016).

For the simulations, the signal of January 2001 was considered in all Monte Carlo samples. The chosen realization of the monthly average of the residual time variable signal (computed as AOHIS–DEAL–AOerr) is shown in figure 16 in terms of both potential and second radial derivative computed at the satellite mean altitude of 360 km.

Zoom In Zoom Out Reset image size

The results obtained in the three simulated scenarios are presented in figures 17 and 18 showing the empirical error degree variance curves for the clock-only and the clock-gradiometer solutions, respectively. Here, the estimation error is compared with the error obtained for the static gravity field only (i.e. following the same simulation scheme used in section 4.2). In order to compute the empirical errors when the signal includes non-tidal temporal variations, the estimates obtained for each Monte Carlo sample are compared with the corresponding reference static model plus the monthly average of the time-variable signal.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The error curves show that there is a degradation in the quality of the solution below degree 20 due to the presence of non-tidal mass variations (compare the solid and dotted lines in figures 17 and 18). As expected, the Bender-class formation significantly improves the capability of detecting the time variable gravity signal with respect to the GRACE-class formation. This is true in the clock-only solution, but especially in the clock-gradiometer joint solution. The transition from two to three in-line satellites per orbit is not really effective in improving the quality of the clock-only solution, while it has a significant impact on the clock-gradiometer joint solution (compare the red and yellow solid lines in figures 17 and 18, respectively). This may be due to the role of the gravity gradients in the solution, but also to the data gridding procedure that is able to locally average the time-variable signal mitigating its impact at a global level.

All the solutions which included non-tidal mass variations were computed by applying the same parameter calibration used in the solutions for the static gravity field (e.g. regularization in the least squares solution, local adapting of covariance function in the gridding procedure, etc, see section 3). Therefore, it can be concluded that part of the error increment shown in figure 18 is related to the space-wise algorithm, which should require a finer tuning to better exploit the time-variable information included in the data.