## 4 Discussion  
According to Figures 1 and 2 as well as Table 4, the derived radial accelerations were noisier than the other two components. This can be explained by a relatively low accuracy of the GPS positioning in the radial direction, which is caused by the geometric configuration of the satellite constellation [44; 45; 46]. Additionally, the along-track accelerations performed weaker than the cross-track component, which was caused by an increase or decrease in the number of visible satellites because geometry changes can inevitably degrade the position solution and create large acceleration errors [46]. Therefore, stringent position accuracies and poor tolerance of geometry changes are the main drawbacks of the orbit differentiation method. However, the carrier phase differentiation method is more robust because no ambiguity resolution is required. In addition, the distribution of gross errors in Figure 1 and geoidal errors in Figures 5 and 7 were geographically correlated and highly inhomogeneous; large errors were mainly located close to the geomagnetic poles. There are two reasons for this situation: the poor observation geometry in geographical polar region, which is very close to the geomagnetic polar region; and the ionospheric scintillation effects, which cause short-period irregular changes in the phase and amplitude of signal [44].  
According to Figures 4 and 6 as well as Tables 5 and 7, the carrier phase solution showed a slightly better performance than the other solutions derived from the orbits only in the high-frequency portion; the 71-day solution corresponds to degree \(n>100\) and 2.5-year solution corresponds to degree \(n>70\)  
Figure 7: Geoid height errors of the recovered spherical harmonic coefficients w.r.t. DGM-1S with the 500 km Gaussian smoothing applied.  
In this frequency band, the acceleration noise of the carrier phase solution was lower than the others, as shown in Figure 2. The degrees 100 and 70 approximately correspond to the frequency 0.019 and 0.013 Hz. In comparison with the Graz-orbit solution, the carrier phase solution had a lower accuracy for the degree \(n<70\) and a higher accuracy for \(n>70\), which is consistent with frequencies less than 0.013 Hz (approximately corresponding to the degree 70). Moreover, the cumulative geoid height differences of the Graz-orbit, ESA-orbit, and carrier phase solutions up to degree 50 were 0.97, 1.01, and 1.08 cm, respectively, as seen in in Table 7, which shows that the Graz-orbit solution had the highest precision for the degrees lower than 50. In Table 5, the carrier phase solution was better than the SGG-orbit solution for almost all degrees (n > 30), and the accuracy improvement of geoid height up to d/o 130 was 21%, which is similar to that in [PERSON] et al. (2015). A comparison of Figures 4 and 6 as well as Tables 5 and 7 indicates that the longer data period has improved the solutions significantly along the entire frequency band. In addition, as depicted in Figures 5 and 7, the pronounced errors along the geomagnetic equator were observed in all solutions. This is a result of the high electron density in addition to large short-term variations in the ionosphere near the geomagnetic equator ([PERSON] et al., 2016; [PERSON] et al., 2016).  
As shown in Figure 3, the zonal and near-zonal spherical harmonics coefficients of the recovered models were worse than the other coefficients due to the fact that the GOCE satellite flies in a sun synchronous orbit with an inclination of 96.7\({}^{\circ}\), which causes the ill-posed problem of the least squares used for recovering the gravity field model ([PERSON] et al., 2016; [PERSON] et al., 2016). However, zonal harmonic coefficients estimated by the carrier phase differentiation method were slightly better than those from the orbit differentiation method.