Impact of Groundtrack Pattern of a Single Pair Mission on the Gravity Recovery Quality

: For future gravity satellite missions, aliasing of high frequency geophysical signals into the lower frequencies is one of the most challenging obstacles to recovering true gravity signals, i.e., to recover the truth. Several studies have investigated the impact of satellite groundtrack pattern on the quality of gravity recovery. Among those works, the concept of sub-cycle has been discussed as well. However, most of that research has focused on the impact of sampling patterns on global solutions up to a ﬁxed maximum spherical harmonic coefﬁcient, rather than the associated coefﬁcient deﬁned by the Colombo-Nyquist and modiﬁed Colombo-Nyquist rules. This work tries to look more closely into the inﬂuence of sampling patterns on the gravity recovery quality for global and regional studies when the Colombo-Nyquist and modiﬁed Colombo-Nyquist rules apply. For the regional study, the impact of groundtrack patterns of different satellite constellation scenarios are investigated for a hydrological basin in central Africa. The quality of the gravity products are assessed by different metrics, e.g., by spatial covariance representation. The potential meaning of the sub-cycle concept in terms of global and local impacts is also investigated by different repeat-orbit scenarios with even and odd parities. Different solution scenarios in terms of the original and modiﬁed Colombo-Nyquist rules will be discussed. The results of our study emphasize the impact of maximum harmonics of the recovery, the inﬂuence of sub-cycle on the local gravity recovery and the mission formation impact on the recovery error.

Furthermore, the European Space Agency (ESA) projects "Assessment of a Next Generation Mission for Monitoring the Variations of Earth's Gravity" [14] and "Assessment of Satellite Constellations for Monitoring the Variations in Earth's Gravity Field" [15] investigated the science requirements, performance criteria and design of future satellite gravity missions. These two ESA studies briefly discussed two advanced formations and Bender configurations. However, because of the complexity of the technical implementation of alternative formations such as Pendulum, Cartwheel and LISA, their employment as a future gravity mission is currently not of interest to the community [11], although the future technical improvements may bring them back to the consideration.
Some of the aforementioned studies also suggested a range of parameter choices for double-pair satellite missions for full repeat period and sub-Nyquist solutions.
Most of the aforementioned works chose individual constellation scenarios by rough assessments of the sampling behavior of the missions. However, the studies by [5,6,8,15] investigated some performance criteria for their optimal double-pair mission scenario search strategies.
In gravity satellite missions, two sampling descriptions exist, which dictate the quality of the gravity recovery change in space and time domains. First, a Heisenberg-like uncertainty theorem in satellite geodesy addresses the trade-off between spatial and temporal resolution by i.e., the product of spatial resolution D space and the time resolution D time is constant and equal to 2πT revolution with T revolution as the revolution time [8,16,17]. Second, the Colombo-Nyquist rule (CNR) requires where β is the number of satellite revolutions within a repeat period in nodal days (α). The formulation means that the number of satellite revolutions in a repeat period should be equal to or larger than twice the maximum spherical harmonic degree (L max ) or order (M max ) to be detected [18,19]. Indeed, this rule limits the spatial resolution of the gravity recovery. However, the works by [20,21] argue that the spatial resolution can be improved, not by twice the maximum degree or order L max or M max , but by a smaller degree or order values. In particular, the work by [21] discusses that those values can even reach L max or M max themselves, depending on the parity of the satellite repeat orbit. However, the description by [21] is not fully applicable, and requires further investigations. In fact, this study states that a recoverable gravity solution in a full repeat period can be achieved by the new modified rule when β − α is odd for the near-polar orbits. Later, [8] took the benefit of the proposed modified Colombo-Nyquist rule (MCNR) to achieve sub-Nyquist gravity recovery, although it should be mentioned that the term "sub-Nyquist" in that work is not defined according to sampling distribution and groundtrack pattern evolution, but by fixed values defined by mission revolutions per nodal days. Therefore, in that context, one can consider Nyquist as the upper limit in sampling theorem. The short-time intervals solutions, called sub-Nyquist solutions, in consequence, may benefit from higher temporal resolution, and thus, less temporal aliasing, compared to the longer span solutions, e.g., defined by CNR. An important benefit of having those solutions is that they can also be employed as de-aliasing products when time-variable gravity solutions of longer time intervals are targeted [22].
Besides the two aforementioned rules, the spatio-temporal groundtrack pattern distribution of a satellite mission scenario can significantly affect the quality of the gravity recovery. Therefore, studying the gap evolution of the satellite mission groundtrack is of a great interest as well. Some studies like [8,[12][13][14] investigate the impact of sampling evolution on the gravity recovery quality. For example, [8,14] look into different categories of orbits where drifting, slow skipping and fast skipping orbits, and their impacts on the quality of gravity solutions for single-pair satellite missions are discussed. Moreover, the work by [12] studies the refinement of Colombo-Nyquist and modified Colombo-Nyquist rules for different latitudes rather than at the equator, where the former studies have focused. The paper by [13], on the other hand, looks into the groundtrack pattern evolution of double pair missions where the influence of relative difference of the right ascension of ascending nodes between the two pairs (∆Ω) is investigated. In fact, the paper discusses three reasons for the quality variations of the gravity solutions: (i) the sampling pattern of each satellite pair of a double pair mission with different repeat periods evolves differently in the time domain.
This choice can therefore address the sampling impact in a clearer way. Furthermore, the particular choice of the RAAN values refers to different sampling behavior, i.e., to the case where sub-cycle happens and to the cases where it does not happen in that area (and at that specific time).

RST Gravity Recovery Simulation Software
For the gravity recovery simulation, a nominal circular orbit simulation package of so-called Reduced-Scale Tool, (RST) (previously called Quick-look Tool, QLT [8,17]) is employed. The software assumes that the satellite orbits are nominal, i.e., the only perturbing force is J 2 effect by the Earth flattening, hence the semi-major axis, inclination angle and eccentricity of the satellites are not time-variable parameters [8]. The observation equation is built up by . ρ 2 (t)) = e 12 .(V(r 2 (t)) − V(r 1 (t))) (3) where ρ is range, . ρ is range rate and ..
ρ is range acceleration. ∆ . r 12 is the relative velocity vector and e 12 is the unit vector of the inter-satellite distance of the satellite pair. The potential gradient at the first and second satellite positions for the epoch t are V(r 1 (t)) and V(r 2 (t)), respectively.
By the RST software, the right side of the equation above is calculated at the nominal positions of the satellites for the time interval of interest. The time-variable potential of the Earth at the positions of the satellites is calculated by the provided time-variable gravity field models at those epochs. Afterwards, the calculated values for the right side of the equation are set to the left side as the observables at those epochs, and the gravitational potential, in terms of spherical harmonics coefficients, is estimated through the system of equations for that time interval. It is important to mention that although the assumption of keeping the satellites in a perfect nominal orbit is not realistic, the simulation software provides a quick comparison of gravity recovery solutions in different mission scenarios, which can later be studied in more details by full-scale recovery tools. An evaluation of the RST has been made with the ll-SST acceleration approach applied to orbits from real orbit integration by [8], which shows a good agreement between the results of the two approaches for gravity recoveries, despite the differences between the two methods.
The equatorial bulge of the Earth (the Earth flattening effect) is responsible for change of the Keplerian elements . ω (argument of perigee rate), . Ω (rate of right ascension of the ascending node) and . M (mean anomaly rate). For circular orbits, the argument of perigee is not defined, but the other two elements change with the following rates: .
where a, e and I are respectively the semi-major axis, eccentricity (here sets to zero for the circular orbits) and inclination angle. The parameter n = GM a 3 stands for the mean motion of the satellite and T is the disturbing potential field.

Geophysical Input Models
For the simulation environment of this study, we employ the dominant mass variations of hydrology, ice and solid Earth of the Earth system by use of the time-variable gravity field models generated within the ESA-project "ESA Earth System Model for Gravity Mission Simulation Studies" [23] as an update to previous ESA models by [24], provided as spherical harmonic coefficients up to degree and order 180 at 6 h intervals. The time-span applied in our study starts from 1 January 1996 for 3-days to 9-days gravity recovery. The difference between two ocean tide models EOT08a [25] and [26], is assumed as the ocean tide error in simulation, where the error is an important source for aliasing of the high frequency signals to the low frequency signals [27]. Moreover, the atmosphere-ocean error (AO error) product from IAPG (TU Munich) is employed for the forward modeling [15]. The AO error model is defined as two atmosphere models difference (ECMWF-NCEP) plus 10% of the ocean signal of the model OMCT. Furthermore, all geophysical forward models are smoothed by Gaussian filter of radius 4000/2L max km to avoid ringing effects (Gibbs phenomenon) caused by the truncation error.

Gap Evolution and the Sub-Cycle Concept
In [14], the sub-cycle concept is defined as follows: Assuming Γ = β/α as the number of orbit revolutions per day, the "Fundamental Interval" S gives the angular space between two Ascending Node Crossing (ANX) consecutive in time The sub-interval S i is the sampling angle after an entire repeat period (RC) Therefore, one can rewrite the parameter S as follows The parameter Γ can also be written as where I is the integer part of Γ. Within S, the ANX of days n and n + 1 are always separated by a distance of N or α − N sub-intervals S i .
Two types of orbits can be recognized: (i) drifting orbits when N = 1 or N = α − 1, and (ii) skipping orbits in the other cases. In a drifting orbit, the sampling of the fundamental interval is very progressive with large unobserved gaps, while a skipping orbit, on the other hand, has a more complex coverage pattern, reducing the persistence of the large unobserved gaps quite fast. In fact, the skipping orbit owns a very wide range of spatial and temporal coverage patterns.
The concept "sub-cycle" (sc) can be defined as the smallest number of days after which an Ascending Node Crossing (ANX) falls at S i or (α − 1)S i from the first ANX. According to this definition, the sub-cycle of a drifting orbit is one day. The sub-cycle can be an interesting parameter to measure how fast an orbit is in reducing the largest unobserved gap [8]. Among different kinds of skipping orbits, some orbits also feature one or more pseudo sub-cycle (psc) which can also be an interesting concept to study the temporal sampling of a satellite scenario [14]. The reader is referred to [8,14] for more details about different categories of orbits and the sub-cycle concept with the gap evolution graphs.

Results
We simulate different inline (GRACE-like) mission scenarios ( Table 1). The scenarios are taken from [8] which try to cover a range of orbits from drifting to fast-skipping with altitudes around 260 to 360 km. The simulations are based on a simplified version of the ll-SST acceleration approach using nominal orbits and the Hill equations [28]. The input models are the relevant time-variable gravity fields of hydrology, ice and solid Earth (HIS) and the difference of two ocean tide models as well as a model for Atmosphere-Ocean (AO) error (see Section 2). The sampling interval of the measurements is ∆t = 1 s. The error is represented as the difference between the gravity recovery (output) and the mean of the HIS time-variable gravity models (input) of the simulations for the same time duration. In all the simulations, the coefficients ∆C 00 , ∆C 10 , ∆C 20 , ∆C 11 and ∆S 11 are neglected, since the ll-SST measurement concepts for these coefficients are poorly sensitive (although depends on the recovery approach). A white noise model of tracking laser interferometer on the level of range acceleration with a PSD (power spectral density) of 10 −10 ms −2 / √ Hz is applied in our RST simulation software.

Global Gravity Recovery Solutions
The global Root Mean Square (RMS) error of the mission scenarios of Table 1 are illustrated in Figures 1-6 for different recovery cases L max = 32, 64, CNR and MCNR. The choices of L max = 32, 64 are arbitrary, but selected as two values below CNR (L max = 32) and between CNR and MCNR (L max = 64) for 6-day solutions (see the next sections). The error is calculated in terms of geoid height for 3-day to 9-day gravity solutions by each mission scenario. As the figures show, the scenario 1 with long full repeat cycle (32 days), but shortest sub-cycle span (1 day), performs the worst. For the MCNR recovery, however, the last scenario (scenario 12) with 8 days repeat cycle, generates extremely large errors on the 8th day. That is clearly not expected, since 8 days is the full repeat cycle of this scenario; therefore, one expects the best performance from the sampling viewpoint for that time-span. The potential reason for such a big error may come from high spatial aliasing of the solution (8-day recovery corresponds to L max = 128 for MCNR) and the geometry of the orbit (short full-repeat cycle scenario).
For the case L max = 32 ( Figure 1), with the increase of the time-pan of the solutions, the recovery error decreases. That is caused by more sampling in space, and therefore, less spatial aliasing, even though the temporal aliasing gradually increases, but not significantly. The same reason also applies for the case L max = 64 (Figures 2 and 3). Figure 2, however, shows extremely high errors for solutions shorter than 4 days. In fact, the error is caused by poor conditioning of the normal matrices (singularity problem). The figure is deliberately brought to show that beyond MCNR (here spherical harmonic, SH, degree 48 for 3-day solutions), the recovery is very poor which is, in particular, the case for some orbit configurations.
For L max = CNR and MCNR, however, it is generally different, i.e., for most of the scenarios, the longer the recovery time-interval, the bigger the error (see . However, different scenarios perform differently. For some scenarios of the case L max = CNR, the error increase is significant. The reason for the increase should be associated with the spatial aliasing (by higher spherical harmonic coefficient degree and order to be solved), which is more predominant for some mission scenarios, most likely caused by the geometry of their orbits, and hence, the time evolution of their groundtracks. A similar statement is also valid for the case MCNR ( Figures 5 and 6), although more fluctuations are seen in the recovery error of few scenarios. In particular, Figure 5 illustrates very large errors for some orbit scenarios. We think that these extreme error values are caused by severe spatial aliasing through  Table 1 mission scenarios for maximum spherical harmonic degree and order 32.    Table 1 mission scenarios for maximum spherical harmonic degree and order 32.  Table 1 mission scenarios for maximum spherical harmonic degree and order 32.             Table 1 mission scenarios for maximum spherical harmonic degree and order MCNR.

Two Scenarios for Comparisons
Here, we choose two mission scenarios from Table 2 for more detailed comparison and investigation. The scenarios have been selected such that they have almost the same altitudes; therefore, the mission height difference does not affect the performance. However, scenario a is in an odd parity orbit, whereas scenario b is in the even parity. Both orbits have almost similar groundtrack gap evolution and same sub-cycle of 6 days (Figure 7), but different full repeat cycles (13 days vs. 31 x   Table 1 mission scenarios for maximum spherical harmonic degree and order MCNR.  Table 1 mission scenarios for maximum spherical harmonic degree and order MCNR.

Two Scenarios for Comparisons
Here, we choose two mission scenarios from Table 2 for more detailed comparison and investigation. The scenarios have been selected such that they have almost the same altitudes; therefore, the mission height difference does not affect the performance. However, scenario a is in an odd parity orbit, whereas scenario b is in the even parity. Both orbits have almost similar groundtrack gap evolution and same sub-cycle of 6 days (Figure 7), but different full repeat cycles (13 days vs. 31 x

Two Scenarios for Comparisons
Here, we choose two mission scenarios from Table 2 for more detailed comparison and investigation. The scenarios have been selected such that they have almost the same altitudes; therefore, the mission height difference does not affect the performance. However, scenario a is in an odd parity orbit, whereas scenario b is in the even parity. Both orbits have almost similar groundtrack gap evolution and same sub-cycle of 6 days (Figure 7), but different full repeat cycles (13 days vs. 31 days). The performance of the mission choices are studied for two different formation flights (GRACE-like and Pendulum) and different angles of right ascension of the ascending node (RAAN). The choices of RAAN for the simulation start date of 1 January 1996 correspond to: -  -Ω = 0°: Starting from a point in Pacific Ocean, west of South-America, -Ω = 116°: Starting from a point in West-Africa, -Ω = 132°: Starting from a point in East-Africa. GRACE-like (Inline) Formation Figure 8 shows the gravity recovery of the GRACE-like scenarios a and b for Lmax = 32.  Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order 32. Different colors stand for different Ω.
As we see in Figure 8, after the 5th day, the gravity recovery by Ω = 132° is better than Ω = 116°. The recovery by Ω = 0° is the worst. The odd parity scenario a usually performs a bit better than the even parity scenario b.
In terms of global recovery, the 6th day as sub-cycle period is not significant. In case L max = 64, Figure 9 shows that the performance of the odd parity scenario is better than the even parity scenario, but the impact of Ω is not significant. In terms of global recovery, the 6th day as sub-cycle period might be considered more significant for L max = 64 than the case of L max = 32.  Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order 64. Different colors stand for different Ω.   Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order 32. Different colors stand for different Ω.
As we see in Figure 8, after the 5th day, the gravity recovery by Ω = 132 • is better than Ω = 116 • . The recovery by Ω = 0 • is the worst. The odd parity scenario a usually performs a bit better than the even parity scenario b.
In terms of global recovery, the 6th day as sub-cycle period is not significant. In case L max = 64, Figure 9 shows that the performance of the odd parity scenario is better than the even parity scenario, but the impact of Ω is not significant.  Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order 32. Different colors stand for different Ω.
As we see in Figure 8, after the 5th day, the gravity recovery by Ω = 132° is better than Ω = 116°. The recovery by Ω = 0° is the worst. The odd parity scenario a usually performs a bit better than the even parity scenario b.
In terms of global recovery, the 6th day as sub-cycle period is not significant. In case L max = 64, Figure 9 shows that the performance of the odd parity scenario is better than the even parity scenario, but the impact of Ω is not significant. In terms of global recovery, the 6th day as sub-cycle period might be considered more significant for L max = 64 than the case of L max = 32.  Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order 64. Different colors stand for different Ω.   Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order 64. Different colors stand for different Ω.
In terms of global recovery, the 6th day as sub-cycle period might be considered more significant for L max = 64 than the case of L max = 32.
For L max = CNR, the dominant effect is the different impacts of odd and even parities on the gravity solutions larger than 7-day ( Figure 10). Moreover, the impact of Ω can be seen (Ω = 132 • usually performs better than Ω = 116 • , and Ω = 0 • is the worst). For Lmax = CNR, the dominant effect is the different impacts of odd and even parities on the gravity solutions larger than 7-day ( Figure 10). Moreover, the impact of Ω can be seen (Ω = 132° usually performs better than Ω = 116°, and Ω = 0° is the worst).
In terms of global recovery, the 6th day as sub-cycle period has a significant influence. Afterwards, the error increases with the increase of time-interval of the solutions, which is very likely caused by the impact of temporal aliasing, since spatial aliasing is tried to be avoided (according to CNR).  Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order CNR. Different colors stand for different Ω.
The case Lmax = MCNR is shown in Figure 11. The figure shows that above 6-day solutions, the odd parity scenario performs better than the even parity. Surprisingly, Ω = 0° performs almost the same as, or even better (for longer time solution) than, the others. Spatial aliasing is thought to be the biggest source of error, which increases as the time period of the gravity solutions increases (although temporal aliasing also exists, as for the other cases). Here, 6-day (as sub-cycle period) does not show important meaning. Possibly, the domination of spatial aliasing and the distribution of the geophysical signals with different strength over the globe hide the meaningful pattern.   Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order CNR. Different colors stand for different Ω.
In terms of global recovery, the 6th day as sub-cycle period has a significant influence. Afterwards, the error increases with the increase of time-interval of the solutions, which is very likely caused by the impact of temporal aliasing, since spatial aliasing is tried to be avoided (according to CNR).
The case L max = MCNR is shown in Figure 11. The figure shows that above 6-day solutions, the odd parity scenario performs better than the even parity. Surprisingly, Ω = 0 • performs almost the same as, or even better (for longer time solution) than, the others. Spatial aliasing is thought to be the biggest source of error, which increases as the time period of the gravity solutions increases (although temporal aliasing also exists, as for the other cases). Here, 6-day (as sub-cycle period) does not show important meaning. Possibly, the domination of spatial aliasing and the distribution of the geophysical signals with different strength over the globe hide the meaningful pattern.

Pendulum Formation
The addition of cross-track measurement of Pendulum formation to only along-track measurement of GRACE-like formation brings the expectation that Pendulum formation should reduce the recovery error through its more isotropic sampling behavior. Figures 12-15 illustrate the gravity recovery of scenarios a and b of Table 2 of a Pendulum formation, respectively for L max = 32, 64, CNR and MCNR. The average intersatellite distance of the two satellites in the formation is 100 km.
For L max = 32 in Figure 12, the recovery by Ω = 132° is better than Ω = 116°. The recovery by Ω = 0° is the worst and gets even worse at the longer periods. No significant difference is seen by the performance of the odd and even parities.
In terms of global recovery, the 6th day as sub-cycle period is not significantly important.    Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order MCNR. Different colors stand for different Ω.

Pendulum Formation
The addition of cross-track measurement of Pendulum formation to only along-track measurement of GRACE-like formation brings the expectation that Pendulum formation should reduce the recovery error through its more isotropic sampling behavior. Figures 12-15 illustrate the gravity recovery of scenarios a and b of Table 2 of a Pendulum formation, respectively for L max = 32, 64, CNR and MCNR. The average intersatellite distance of the two satellites in the formation is 100 km.
For L max = 32 in Figure 12, the recovery by Ω = 132 • is better than Ω = 116 • . The recovery by Ω = 0 • is the worst and gets even worse at the longer periods. No significant difference is seen by the performance of the odd and even parities. In terms of global recovery, the 6th day as sub-cycle period is not significantly important.  Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order MCNR. Different colors stand for different Ω.

Pendulum Formation
The addition of cross-track measurement of Pendulum formation to only along-track measurement of GRACE-like formation brings the expectation that Pendulum formation should reduce the recovery error through its more isotropic sampling behavior. Figures 12-15 illustrate the gravity recovery of scenarios a and b of Table 2 of a Pendulum formation, respectively for L max = 32, 64, CNR and MCNR. The average intersatellite distance of the two satellites in the formation is 100 km.
For L max = 32 in Figure 12, the recovery by Ω = 132° is better than Ω = 116°. The recovery by Ω = 0° is the worst and gets even worse at the longer periods. No significant difference is seen by the performance of the odd and even parities. In terms of global recovery, the 6th day as sub-cycle period is not significantly important.    Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order 32. Different colors stand for different Ω. Figure 13 illustrates better performance in recovery quality by Ω = 132 • than Ω = 116 • . After 5-day, the gravity recovery by Ω = 0 • is the worst, and gets even worse for the longer periods. No significant difference is seen by the performance of the odd and even parities, although the even parity performs only a bit better than the odd parity after 5-day. In terms of global recovery, the 6th day as sub-cycle period is not significantly important, although it shows a small impact on the recovery by Ω = 0 • and Ω = 116 • .  Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order 32. Different colors stand for different Ω. Figure 13 illustrates better performance in recovery quality by Ω = 132° than Ω = 116°. After 5day, the gravity recovery by Ω = 0° is the worst, and gets even worse for the longer periods. No significant difference is seen by the performance of the odd and even parities, although the even parity performs only a bit better than the odd parity after 5-day. In terms of global recovery, the 6th day as sub-cycle period is not significantly important, although it shows a small impact on the recovery by Ω = 0° and Ω = 116°.  Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order 64. Different colors stand for different Ω.
Unlike the case of CNR by the GRACE-like formation, the errors for Ω =116° and Ω = 132° in Figure 14 decrease by the longer period of the solutions. One may associate this error reduction to the less spatial aliasing. However, that is not the case for Ω = 0°. The question is then how the combination of isotropy and sampling plays role.
Looking at the Figure 14, the Ω = 132° scenario shows better performance than Ω = 116° which outperforms Ω = 0°. Furthermore, the even parity performs a bit better than the odd parity.
In terms of global recovery, the 6th day as sub-cycle period is not significantly important.   Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order 64. Different colors stand for different Ω.
Unlike the case of CNR by the GRACE-like formation, the errors for Ω = 116 • and Ω = 132 • in Figure 14 decrease by the longer period of the solutions. One may associate this error reduction to the less spatial aliasing. However, that is not the case for Ω = 0 • . The question is then how the combination of isotropy and sampling plays role.
Looking at the Figure 14, the Ω = 132 • scenario shows better performance than Ω = 116 • which outperforms Ω = 0 • . Furthermore, the even parity performs a bit better than the odd parity.
In terms of global recovery, the 6th day as sub-cycle period is not significantly important.  In Figure 15 (MCNR case), the error increases for all cases after 6-day solution. In this case, extracting a meaning for sub-cycle (6-day) is not straightforward.
The odd parity outperforms the even parity at longer periods of solutions. The case Ω = 132° usually performs better than the two other cases.  Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order MCNR. Different colors stand for different Ω.

Covariance Matrices
The illustration of the covariance matrices of the gravity solutions by different scenarios and cases can be employed as an indicator for the associated error levels, very much influenced by the geometry of the orbits and their sampling behaviors.  Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order CNR. Different colors stand for different Ω.
In Figure 15 (MCNR case), the error increases for all cases after 6-day solution. In this case, extracting a meaning for sub-cycle (6-day) is not straightforward.
The odd parity outperforms the even parity at longer periods of solutions. The case Ω = 132 • usually performs better than the two other cases.  In Figure 15 (MCNR case), the error increases for all cases after 6-day solution. In this case, extracting a meaning for sub-cycle (6-day) is not straightforward.
The odd parity outperforms the even parity at longer periods of solutions. The case Ω = 132° usually performs better than the two other cases.  Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order MCNR. Different colors stand for different Ω.

Covariance Matrices
The illustration of the covariance matrices of the gravity solutions by different scenarios and cases can be employed as an indicator for the associated error levels, very much influenced by the geometry of the orbits and their sampling behaviors.   Table 2, respectively) for the gravity solutions up to maximum spherical harmonic degree and order MCNR. Different colors stand for different Ω.

Covariance Matrices
The illustration of the covariance matrices of the gravity solutions by different scenarios and cases can be employed as an indicator for the associated error levels, very much influenced by the geometry of the orbits and their sampling behaviors. Figures 16 and 17 depict the covariance matrices of 6-day solutions of the scenario (a) inline and Pendulum formations for Ω = 116 • , up to L max = 32, L max = 64, L max = CNR and L max = MCNR solutions. The illustrations imply that for the inline formation when the maximum spherical harmonic degree and order goes above the CNR rule for 6-day (L max = 64 and L max = MCNR), the level of error shown by the arcs and off-diagonal elements (due to the spatial aliasing) gets significant. For the Pendulum formation, however, those error patterns are not as strong as the inline formation. The illustrations imply that for the inline formation when the maximum spherical harmonic degree and order goes above the CNR rule for 6-day (Lmax = 64 and Lmax = MCNR), the level of error shown by the arcs and off-diagonal elements (due to the spatial aliasing) gets significant. For the Pendulum formation, however, those error patterns are not as strong as the inline formation.

Inline (GRACE-like) Formation
Inline (GRACE-like) Formation The illustrations imply that for the inline formation when the maximum spherical harmonic degree and order goes above the CNR rule for 6-day (Lmax = 64 and Lmax = MCNR), the level of error shown by the arcs and off-diagonal elements (due to the spatial aliasing) gets significant. For the Pendulum formation, however, those error patterns are not as strong as the inline formation.

West Africa as an example
Two scenarios of Table 2 with almost same altitude and similar gap evolution (both has 6-day sub-cycle time-span) are investigated for their regional gravity recovery. Figure 18 depicts the groundtrack of the two scenarios on the selected area in West Africa for different RAAN angles after 6 days.

West Africa as an Example
Two scenarios of Table 2 with almost same altitude and similar gap evolution (both has 6-day sub-cycle time-span) are investigated for their regional gravity recovery. Figure 18 depicts the groundtrack of the two scenarios on the selected area in West Africa for different RAAN angles after 6 days.

West Africa as an example
Two scenarios of Table 2 with almost same altitude and similar gap evolution (both has 6-day sub-cycle time-span) are investigated for their regional gravity recovery. Figure 18 depicts the groundtrack of the two scenarios on the selected area in West Africa for different RAAN angles after 6 days.  The error RMS of the gravity recovery in terms of geoid height for global and regional (West-Africa) 6-day solution of different cases are shown in Tables 3-10. The results from the tables imply that the best global performance (the smallest error) usually happens at Ω = 132°. That might be expected, since among the three different Ω values, the mission track passes over the biggest signals when Ω = 132° (i.e., over East-Africa). In terms of regional investigation, the results are usually optimal when Ω = 116°. That is, in fact, the area where sub-cycle happens. However, for both global and regional cases, the above statements are not always valid. When the maximum spherical degree and order of the solution goes above the Colombo-Nyquist rule (CNR), exceptions can be observed. One possible explanation for those exceptions is the large spatial aliasing and the leakage from the surroundings which may deteriorates the impact of sub-cycle. Clearly, more detailed studies should be conducted in this direction.
The impact of the mission formation on the error can also be seen in the results. The outcomes show that the Pendulum formation usually outperforms the inline (GRACE-like) formation. Table 3. Error RMS in terms of geoid height for global and regional (West-Africa) 6-day recovery up to Lmax = 32 of the scenario (a) inline and Pendulum formations. The smallest values in global and regional cases for each inline and Pendulum formation are written in blue.

Inline (GRACE-Like)
Pendulum  Table 4. Error RMS in terms of geoid height for global and regional (West-Africa) 6-day recovery up to Lmax = 32 of the scenario (b) inline and Pendulum formation. The smallest values in global and regional cases for each inline and Pendulum formation are written in blue.

Inline (GRACE-Like)
Pendulum The error RMS of the gravity recovery in terms of geoid height for global and regional (West-Africa) 6-day solution of different cases are shown in Tables 3-10. The results from the tables imply that the best global performance (the smallest error) usually happens at Ω = 132 • . That might be expected, since among the three different Ω values, the mission track passes over the biggest signals when Ω = 132 • (i.e., over East-Africa). In terms of regional investigation, the results are usually optimal when Ω = 116 • . That is, in fact, the area where sub-cycle happens. However, for both global and regional cases, the above statements are not always valid. When the maximum spherical degree and order of the solution goes above the Colombo-Nyquist rule (CNR), exceptions can be observed. One possible explanation for those exceptions is the large spatial aliasing and the leakage from the surroundings which may deteriorates the impact of sub-cycle. Clearly, more detailed studies should be conducted in this direction.
The impact of the mission formation on the error can also be seen in the results. The outcomes show that the Pendulum formation usually outperforms the inline (GRACE-like) formation. Table 3. Error RMS in terms of geoid height for global and regional (West-Africa) 6-day recovery up to L max = 32 of the scenario (a) inline and Pendulum formations. The smallest values in global and regional cases for each inline and Pendulum formation are written in blue.

Inline (GRACE-Like)
Pendulum  Table 4. Error RMS in terms of geoid height for global and regional (West-Africa) 6-day recovery up to L max = 32 of the scenario (b) inline and Pendulum formation. The smallest values in global and regional cases for each inline and Pendulum formation are written in blue.  Table 5. Error RMS in terms of geoid height for global and regional (West-Africa) 6-day recovery up to L max = 64 of the scenario (a) inline and Pendulum formations. The smallest values in global and regional cases for each inline and Pendulum formation are written in blue.  Table 6. Error RMS in terms of geoid height for global and regional (West-Africa) 6-day recovery up to L max = 64 of the scenario (b) inline and Pendulum formations. The smallest values in global and regional cases for each inline and Pendulum formation are written in blue.  Table 7. Error RMS in terms of geoid height for global and regional (West-Africa) 6-day recovery up to L max = CNR of the scenario (a) inline and Pendulum formations. The smallest values in global and regional cases for each inline and Pendulum formation are written in blue.  Table 8. Error RMS in terms of geoid height for global and regional (West-Africa) 6-day recovery up to L max = CNR of the scenario (b) inline and Pendulum formations. The smallest values in global and regional cases for each inline and Pendulum formation are written in blue.  Table 9. Error RMS in terms of geoid height for global and regional (West-Africa) 6-day recovery up to L max = MCNR of the scenario (a) inline and Pendulum formations. The smallest values in global and regional cases for each inline and Pendulum formation are written in blue.  Table 10. Error RMS in terms of geoid height for global and regional (West-Africa) 6-day recovery up to L max = MCNR of the scenario (b) inline and Pendulum formations. The smallest values in global and regional cases for each inline and Pendulum formation are written in blue.

Spatial Covariance
Another interesting tool to investigate the sampling impact is spatial covariance illustrations [29]. One may expect that denser sampling should result in more symmetric spatial covariance pattern.
The spatial covariance patterns, centered on the selected points of Figures 19 and 20 (marked by green star at longitude 17 • and latitude 2.5 • ) for inline and Pendulum formations are depicted in Figures 21-28 where two different RAAN angles (Ω = 0 • and Ω = 116 • ) and different maximum spherical harmonic coefficient degrees and orders are addressed for 6-day solutions. The calculation of the spatial covariance is based on the work by [29].
The depicted results of the inline formation show the symmetry of the spatial covariance around the selected points for the cases L max = 32, 64 and CNR, although in the case of L max = 64 the spatial covariance graph center is shifted a bit from the center (in this case, the Ω = 116 • , where sub-cycle happens, shows less affected than Ω = 0 • ). For the Pendulum formation, the graphs of L max = 32 and CNR show symmetric feature, but the graph centers are shifted towards left (more predominant for L max = 32). That is not the case for L max = 64. Interestingly, for the MCNR, the scenario a significantly outperforms the scenario b in terms of the symmetry.  Table 10. Error RMS in terms of geoid height for global and regional (West-Africa) 6-day recovery up to Lmax = MCNR of the scenario (b) inline and Pendulum formations. The smallest values in global and regional cases for each inline and Pendulum formation are written in blue.

Spatial Covariance
Another interesting tool to investigate the sampling impact is spatial covariance illustrations [29]. One may expect that denser sampling should result in more symmetric spatial covariance pattern.
The spatial covariance patterns, centered on the selected points of Figures 19 and 20 (marked by green star at longitude 17° and latitude 2.5°) for inline and Pendulum formations are depicted in Figures 21-28 where two different RAAN angles (Ω = 0° and Ω = 116°) and different maximum spherical harmonic coefficient degrees and orders are addressed for 6-day solutions. The calculation of the spatial covariance is based on the work by [29].
The depicted results of the inline formation show the symmetry of the spatial covariance around the selected points for the cases Lmax = 32, 64 and CNR, although in the case of Lmax = 64 the spatial covariance graph center is shifted a bit from the center (in this case, the Ω = 116°, where sub-cycle happens, shows less affected than Ω = 0°). For the Pendulum formation, the graphs of Lmax = 32 and CNR show symmetric feature, but the graph centers are shifted towards left (more predominant for Lmax = 32). That is not the case for Lmax = 64. Interestingly, for the MCNR, the scenario a significantly outperforms the scenario b in terms of the symmetry.

Discussion
This work aimed to look more closely into the sampling patterns impact on the gravity recovery quality for global and regional studies, as the Colombo-Nyquist and modified Colombo-Nyquist rules apply. We tried to focus on the pure sampling impact in this paper, so no post-processing was considered. In reality, smoothing global or local filter (as one of the filter choices) is implemented. However, higher quality raw solutions (for example, due to denser sampling, and hence less spatial aliasing) can result in using more relaxed filters (e.g., a Gaussian smoothing with smaller radius). That would be then very beneficial in terms of signal-to-noise ratio, although the subject is outside the scope of this paper.
In particular, for the regional study, the groundtrack patterns impact of different satellite constellation scenarios have been investigated for a hydrological basin in central Africa. The quality of the gravity products have been assessed by different metrics such as spatial covariance representation. We also tried to investigate the potential meaning of the sub-cycle concept in terms of global and local impacts by different repeat-orbit scenarios with even and odd parities. Furthermore, different recovery scenarios in terms of the original and modified Colombo-Nyquist rules have been discussed.
The main points of our results can be summarized as follows: In general, for 12 selected mission scenarios of Table 1: For the case L max = 32, with the increase of the time-pan of the solutions, the recovery error decreases. This is because of more sampling in space, and hence, less spatial aliasing, although the temporal aliasing gradually increases (but not significantly). The same argument applies for the case L max = 64 as well. For L max = CNR and MCNR, it is generally different. That means that for most of the scenarios, the longer the recovery time-span, the bigger the error. But here, it is important to notice that different scenarios have different performances. For some scenarios of the case L max = CNR, the error increase is significant. The cause of such an increase should be associated with the spatial aliasing which is more predominant for some mission scenarios, most likely because of the geometry of their orbits and therefore their time evolution of their groundtracks. A similar statement also applies for the case MCNR, even though more fluctuations are seen in the recovery error of few scenarios.
In terms of global recovery and for the selected inline formation scenarios of Table 2: Considering the inline formation missions, for L max = 32, the recovery by Ω = 132 • after the 5th day is better than Ω = 116 • . The recovery by Ω = 0 • is the worst. Odd parity usually performs a bit better than the even parity. The 6th day as sub-cycle period is not significantly important in the global point of view. However, the concept might be more noticeable for L max = 64. For this case, the performance of the odd parity is better than the even parity. That is also the case for L max = CNR. In this case, different impact of odd and even parities on the gravity solutions larger than 7-day are significant. Moreover, the impact of Ω can be seen (Ω = 132 • usually performs better than Ω = 116 • , and Ω = 0 • is the worst). In terms of global recovery, the 6th day as sub-cycle period also has a significant influence. After the 6th day, the error increases with the increase of the recovery time-interval. This effect is mainly assumed to be associated with the impact of temporal aliasing, since the spatial aliasing is tried to be largely avoided, according to CNR (although the spatial aliasing still exits by the impact of orbit geometry).
For the MCNR recovery, the odd parity performs better than the even parity above 6-day solutions. Surprisingly, Ω = 0 • performs almost equal or even better (for longer time recovery) than the other two cases. Spatial aliasing is thought to be the biggest source of the error which increases as time span of the solutions increases (although temporal aliasing also contributes, as it is valid for all the other cases).
6-day (as sub-cycle period) does not show noticeable meaning. Very likely, the domination of spatial aliasing and the distribution of the geophysical signals with different strength over the globe hinder the emergence of any meaningful pattern.
Considering the Pendulum formation of the Table 2 scenarios, and in the terms of global recovery, for the L max = 32, Ω = 132 • performs better than Ω = 116 • . The recovery by Ω = 0 • is the worst and getting even worse over longer periods. One can not see significant difference by the performance of the odd and even parities. The 6th day as sub-cycle period is not significantly important in this case. For L max = 64, the recovery by Ω = 132 • is better than Ω = 116 • . After day 5, recovery by Ω = 0 • is the worst and getting even worse at the longer periods. No significant difference is seen by the performance of the odd and even parities, although the even parity performs only a bit better than the odd parity after day 5. In terms of global recovery, the 6th day as sub-cycle period shows no important feature, although it may depict a very small impact on the recovery by Ω = 0 • and Ω = 116 • .
Unlike the case of CNR by inline formation, reduction of errors for Ω =116 • and Ω = 132 • by the longer period of the solutions are seen (possibly by less spatial aliasing). However, that is not the case for Ω = 0 • . Therefore, one potential question to be addressed in future works might be the way that combination of isotropy (by Pendulum formation) and sampling plays role. The Ω = 132 • case performs better than Ω = 116 • which outperforms Ω = 0 • . The even parity illustrates a slightly better performance than the odd parity. In terms of global recovery, the 6th day (as sub-cycle period) should not be considered significantly important.
For the MCNR case, the error increases for all cases after the 6-day solution. Therefore, it is not clear whether day 6 has meaning in terms of being sub-cycle or not. It is noticeable to see the better performance of the odd parity than the even parity at longer periods of solutions. The case Ω = 132 • mostly performs better than the two other cases.
This paper also illustrates covariance matrices for 6-day solutions of the scenario (a) inline and Pendulum formations for Ω = 116 • , where L max = 32, L max = 64, L max = CNR and L max = MCNR. The illustrations show that for the inline formation when the maximum spherical harmonic coefficient degree and order are above the CNR rule for 6-day (i.e., L max = 64 and L max = MCNR), the error level gets significant (assumed to be associated with spatial aliasing). However, for the Pendulum formation, those error patterns are not as noticeable as the inline formation (very likely associated to the more isotropic sampling by Pendulum formation).
Interpreting the results above, it should be emphasized that with our starting date choice for recovery simulation (1 January 1996), Ω = 0 • corresponds to the groundtrack starting point at Pacific Ocean, east of South America, while Ω = 116 • and Ω = 132 • respectively stand for starting points at West-Africa and East-Africa. The results of this paper show that the best global performance (the smallest error) usually happens at Ω = 132 • . This outcome might be expected since, among the three different values of Ω, the satellite mission track passes over significant hydrological signals when Ω = 132 • (i.e., over East-Africa and the corresponding regions in the same longitudes). In terms of regional investigation, the results are usually optimal for the Ω = 116 • case. In fact, this outcome is also expected, since that is the area where sub-cycle happens for the region targeted in this study (West-Africa). Nevertheless, for both global and regional study, the above statements are not always valid. When the maximum spherical degree and order of the solution goes beyond the Colombo-Nyquist rule (CNR), exceptions can be observed. A possible explanation for those exceptions could be the large spatial aliasing and the leakage from the surroundings which very likely deteriorates the impact of the sub-cycle. Obviously, more detailed studies should be directed in this way.
The results of our study also emphasize the impact of the mission formation on the error, where the Pendulum formation usually outperforms the inline (GRACE-like) formation, as expected through its more isotropic sampling behavior.
The spatial covariance illustrations of this research work for the GRACE-like formation depict the symmetry of the spatial covariance around the selected points for the cases L max = 32, 64, and CNR, although in the case of L max = 64, the spatial covariance graph center is shifted a bit, although for the Ω = 116 • (where sub-cycle happens) the shift is less noticeable than for Ω = 0 • . For the Pendulum formation, the graphs of L max = 32 and CNR feature symmetric patterns, but the graph centers are shifted towards left where it is more predominant for L max = 32. That is not the case for maximum spherical harmonic coefficient 64. It is interesting to say that for the MCNR recovery, the scenario (a) significantly outperforms the scenario (b) in terms of the symmetry.
Regarding the orbit parity impact on the gravity recovery, one should notice that the parity state (odd or even) is usually meaningful for the full repeat cycle. However, the evolution of the groundtrack pattern before the full cycle might be sophisticated and differ in each case; hence it is more reasonable to see the impact by its association to the groundtrack pattern distribution at a specific time span. That means, for the recovery time-spans shorter than the full-repeat cycle, instead of talking about the parity impact on the gravity recovery, it makes more sense to study the groundtrack pattern at that specific time interval and look for its potential association with the recovery quality. This issue will be addressed in our future works in more detail.

Conclusions
This paper examines the impact of sampling patterns on gravity recovery quality for global and regional studies when the Colombo-Nyquist and modified Colombo-Nyquist rules apply. For the regional study, the groundtrack pattern impact of different satellite constellation scenarios have been studied for a hydrological basin in central Africa. The quality of the gravity products were assessed by different metrics, e.g., by spatial covariance representation. The possible meaning of the sub-cycle concept in terms of global and local impacts was also investigated by different repeat-orbit scenarios with even and odd parities where different solution scenarios in terms of the original and modified Colombo-Nyquist rules have been discussed. Our results emphasize the impact of maximum harmonics of the recovery (both global and regional), the influence of sub-cycle on local gravity recovery, and the mission formation impact on recovery error.