# Combination Analysis of Future Polar-Type Gravity Mission and GRACE Follow-On
## 2 Numerical Simulation Scheme
### Mathematical Model for Gravity Field Recovery
Substituting \(\mathbf{r}(t)=\mathbf{r}^{OBS}(t)+\mathbf{v}_{r}(t)\) into the first equation of Equation (4), we get the observational equation for position at epoch \(t\),  
\[\mathbf{v}_{r}(t)=\frac{\partial\mathbf{r}^{REF}(t)}{\partial\mathbf{x}_{0}}\delta\mathbf{x}- \left[\mathbf{r}^{OBS}(t)-\mathbf{r}^{REF}(t)\right]=\frac{\partial\mathbf{r}^{REF}(t)}{ \partial\mathbf{x}_{0}}\delta\mathbf{x}-\Delta\mathbf{r}(t) \tag{5}\]  
where \(\mathbf{r}^{OBS}(t)\) and \(\mathbf{v}_{r}(t)\) represent the position observation and its correction at epoch \(t\). The inter-satellite range rate \(\dot{\rho}(t)\) between GRACE-A and B can be expressed as,  
\[\dot{\rho}(t)=\mathbf{e}_{AB}^{T}(t)\dot{\mathbf{r}}_{AB}(t) \tag{6}\]  
where, the subscripts A and B denote the values of GRACE-A and B, \(\dot{\mathbf{r}}_{AB}(t)=\dot{\mathbf{r}}_{B}(t)-\dot{\mathbf{r}}_{A}(t)\) is the velocity difference vector and \(\mathbf{e}_{AB}(t)=\mathbf{r}_{AB}(t)/\sqrt{\mathbf{r}_{AB}^{T}(t)\mathbf{r}_{AB}(t)}\) with \(\mathbf{r}_{AB}=\mathbf{r}_{B}-\mathbf{r}_{A}\). By substituting \(\dot{\rho}(t)=\dot{\rho}^{OBS}(t)+v_{\dot{\rho}}(t)\) and the second equation of Equation (4) into Equation (6) for both GRACE-A and B, we get the following observational equation for range rate measurement at epoch \(t\),  
\[v_{\dot{\rho}}(t)=\frac{\partial\dot{\rho}^{REF}(t)}{\partial\mathbf{x}_{0}}\delta \mathbf{x}-\left[\dot{\rho}^{OBS}(t)-\dot{\rho}^{REF}(t)\right]=\frac{\partial\dot {\rho}^{REF}(t)}{\partial\mathbf{x}_{0}}\delta\mathbf{x}-\Delta\dot{\rho}(t) \tag{7}\]  
where, \(\dot{\rho}^{REF}(t)\) is computed with the reference orbit, \(\dot{\rho}^{OBS}(t)\) and \(v_{\dot{\rho}}(t)\) stand for the range rate observation and its correction, respectively. The partial derivatives \(\frac{\partial\mathbf{r}^{REF}(t)}{\partial\mathbf{x}_{0}}\) and \(\frac{\partial\mathbf{r}^{REF}(t)}{\partial\mathbf{x}_{0}}\) are numerically integrated by solving the so-called 'variational equation' and \(\frac{\partial\dot{\rho}^{REF}(t)}{\partial\mathbf{x}_{0}}\) is the combination of them, which can be referred to [50] (pp. 240-243), [51] (pp. 30-32) and [49] for detailed descriptions. The observational equations in the form of Equations (5) and (7) for all epochs within one arc, i.e., arc \(j\), are collected together and briefly expressed in the following form,  
\[\mathbf{v}_{j}=\mathbf{A}_{j}\delta\mathbf{x}_{j}-\mathbf{y}_{j}=\left( \begin{array}{cc}\bm {A}_{\mathbf{w}_{j}}&\mathbf{A}_{\mathbf{u}_{j}}\end{array} \right)\left( \begin{array}{ cc}\delta\mathbf{w}_{j}\\ \delta\mathbf{u}\end{array} \right)-\mathbf{y}_{j} \tag{8}\]  
where \(\mathbf{v}_{j}\) denotes the vector of all corrections, \(\mathbf{A}_{j}\), \(\mathbf{A}_{\mathbf{w}_{j}}\) and \(\mathbf{A}_{\mathbf{u}_{j}}\) are the design matrices, and \(\mathbf{y}_{j}\) is the vector of all constant terms of the arc \(j\). The Normal Equation (NEQ) of the arc based least squares adjustment is derived as [48],  
\[\left( \begin{array}{cc}\mathbf{N}_{\mathbf{w}_{j}\mathbf{w}_{j}}&\mathbf{N}_{\mathbf{w}_{j}\mathbf{ u}_{j}}\\ \mathbf{N}^{T}_{\mathbf{w}_{j}\mathbf{u}_{j}}&\mathbf{N}_{\mathbf{u}_{j}\mathbf{u}_{j}}\end{array}  \right)\left( \begin{array}{cc}\delta\mathbf{w}_{j}\\ \delta\mathbf{u}\end{array} \right)=\left( \begin{array}{cc}\mathbf{R}_{\mathbf{w}_{j}} \\ \mathbf{R}_{\mathbf{u}_{j}}\end{array} \right) \tag{9}\]  
where \(\mathbf{N}_{\mathbf{w}_{j}\mathbf{w}_{j}}=\mathbf{A}^{T}_{\mathbf{w}_{j}}\mathbf{P}_{j}\mathbf{A}_\mathbf{{w}_{j}}\), \(\mathbf{N}_{\mathbf{u}_{j}\mathbf{u}_{j}}=\mathbf{A}^{T}_{\mathbf{u}_{j}}\mathbf{P}_{j}\mathbf{A}_\mathbf{{u}_{j}}\), \(\mathbf{N}_{\mathbf{w}_{j}\mathbf{u}_{j}}=\mathbf{A}^{T}_{\mathbf{w}_{j}}\mathbf{P}_{j}\mathbf{A}_\mathbf{{u}_{j}}\) and \(\mathbf{R}_{\mathbf{w}_{j}}=\mathbf{A}^{T}_{\mathbf{w}_{j}}\mathbf{P}_{j}\mathbf{y}_{j}\), \(\mathbf{R}_{\mathbf{u}_{j}}=\mathbf{A}^{T}_{\mathbf{u}_{j}}\mathbf{P}_{j}\mathbf{y}_{j}\), \(\mathbf{P}_{j}\) denotes the weight matrix constructed based on the precision of observations.