## 2 Methods
### Point-Wise Acceleration Approach for Recovering Gravity Field Model  
The point-wise acceleration approach for recovering the gravity field model is based on [PERSON]'s equation of motion in the inertial reference frame [18]:  
\[\ddot{\mathbf{r}}-\mathbf{a}_{cf}-\mathbf{a}_{ncf}=R_{E}^{I}\
abla V(r,\theta,\lambda) \tag{13}\]  
where \(\ddot{\mathbf{r}}\) is the satellite acceleration; \(\mathbf{a}_{cf}\) denotes the time-variable gravity field signals; \(\mathbf{a}_{ncf}\) is the non-gravitational acceleration; \(R_{E}^{I}\) denotes the rotation matrix from the Earth-fixed reference frame to the inertial reference frame; \(\
abla\) is regarded as the gradient symbol; and \(V\) is the static gravitational potential.  
Although the GOCE satellite makes use of a drag-free control system in the flight direction that has largely compensated for the non-gravitational perturbation (mainly air-drag), it is indispensable to accommodate for the systematic errors caused by satellite attitude control, the resulting movement of center-of-mass, and the inconsistent derivation scale factors of each accelerometers. 1-cpr empirical acceleration parameters were often used to absorb the remaining non-gravitational accelerations [30]. We also estimated SH coefficients by estimating the 1 cpr empirical parameters. There was not a significant difference in the solution without estimating the empirical parameters. Accordingly, we preferred to determine SH coefficients without estimating the empirical parameters. All time-variable gravity field signals that include third-body accelerations, tidal accelerations, and relativistic effects were computed according to background models [31; 32; 33; 34].  
\begin{table}
\begin{tabular}{c c c c c} \hline \hline
**Component** & **Min** & **Max** & **Mean** & **RMS** \\ \hline X & \(-3.968\times 10^{-5}\) & \(3.809\times 10^{-5}\) & \(4.522\times 10^{-10}\) & \(3.039\times 10^{-7}\) \\ Y & \(-1.612\times 10^{-5}\) & \(1.563\times 10^{-5}\) & \(-1.649\times 10^{-10}\) & \(1.105\times 10^{-7}\) \\ Z & \(-1.457\times 10^{-5}\) & \(1.313\times 10^{-5}\) & \(1.978\timesAs described in [PERSON] and [PERSON] [35], the gravitational potential \(V\) is usually represented in spherical harmonic series expansion:  
\[V(r,\theta,\lambda)=\frac{GM}{R}\sum_{n=0}^{N}\left(\frac{R}{r}\right)^{n+1}\sum _{m=0}^{n}\left(\overline{C}_{nm}\cos m\lambda+\overline{S}_{nm}\sin m\lambda \right)\overline{P}_{nm}(\cos\theta) \tag{14}\]  
where \((r,\theta,\lambda)\) denotes the spherical polar coordinates with radius, co-latitude, and longitude; \(GM\) is the geocentric gravitational constant; \(R\) is the Earth's mean radius; \(\overline{P}_{nm}\) denotes the fully normalized associated Legendre functions of degree \(n\) and order \(m\); \(N\) is the maximum spherical harmonic degree; and \(\overline{C}_{nm}\) and \(\overline{S}_{nm}\) are the geopotential coefficients. Combining Equations (13) and (14), a regression model can be formed:  
\[\mathbf{y}=B\boldsymbol{x}+\Delta \tag{15}\]  
Therein, \(\mathbf{y}\) denotes the observation vectors; \(\Delta\) denotes the error vectors; \(B\) is the design matrix; and \(\boldsymbol{x}\) is the unknown geopotential coefficient.  
Using Equation (15), the geopotential coefficients can be estimated based on the least-squares method from the satellite accelerations. To guarantee an unbiased estimate, we did not use a regularization technique in the least squares.