A Comparative Study on the Solar Radiation Pressure Modeling in GPS Precise Orbit Determination

For Global Positioning System (GPS) precise orbit determination (POD), the solar radiation pressure (SRP) is the dominant nongravitational perturbation force. Among the current SRP models, the ECOM and box-wing models are widely used in the International GNSS Service (IGS) community. However, the performance of different models varies over different GPS satellites. In this study, we investigate the performances of different SRP models, including the box-wing and adjustable box-wing as a priori models, and ECOM1 and ECOM2 as parameterization models, in the GPS POD solution from 2017 to 2019. Moreover, we pay special attention to the handling of the shadow factor in the SRP modeling for eclipsing satellites, which is critical to achieve high-precision POD solutions but has not yet been fully investigated. We demonstrate that, as an a priori SRP model, the adjustable box-wing has better performance than the box-wing model by up to 5 mm in the orbit day boundary discontinuity (DBD) statistics, with the largest improvement observed on the BLOCK IIR satellites using the ECOM1 as a parameterization SRP model. The box-wing model shows an insignificant orbit improvement serving as the a priori SRP model. For the eclipsing satellites, the three-dimensional (3D) root mean square (RMS) values of orbit DBD are improved when the shadow factor is applied only in the D direction (pointing toward to Sun) than that in the three directions (D, Y, and B) in the satellite frame. Different SRP models have comparable performance in terms of the Earth rotation parameter (ERP) agreement with the IERS EOP 14C04 product, whereas the magnitude of the length of day (LoD) annual signal is reduced when the shadow factor is applied in the D direction than in the three directions. This study clarifies how the shadow factor should be applied in the GPS POD solution and demonstrates that the a priori adjustable box-wing model combined with ECOM1 is more suitable for high-precision GPS POD solutions, which is useful for the further GNSS data analysis.


Introduction
The Global Positioning System (GPS) serves as a critical role in the positioning, navigation, and timing (PNT) services and is also a powerful tool in monitoring and understanding the geophysical processes of the Earth system [1,2], crust movement [3,4], and meteorological studies [5,6]. The International GNSS Service (IGS) community has been making continuous efforts in understanding the errors and further improving satellite orbit dynamic modeling. Currently, the GPS orbit products from IGS and other analysis lefts (ACs) have a one-dimensional (1D) mean error of around 1 cm in terms of inter-AC agreement [7] and 2.5 cm root mean square (RMS) in terms of the satellite laser ranging (SLR) validation [http://www.igs.org (accessed on 9 January 2021)]. One of the major nongravitational forces affecting the current GPS orbits is the solar radiation pressure (SRP), which relies on the knowledge of satellite-related properties. Although many SRP models have been developed, the performance of different SRP models have not been fully investigated, especially considering the different type of GPS satellites and the performances in eclipse seasons [8]. Note that the eclipse season in this study means that the β angle, that is, the elevation angle of the Sun over orbital plane, is within ±13.25 • [9].
The SRP models are mainly divided into two categories, the analytical and empirical models. Analytical models such as the box-wing and ray tracing approach are established based on the physical properties of the satellite components, including the area, reflection, diffusion, absorption, and the satellite attitude. The attitude control modes of GPS satellites are known clearly except for BLOCK IIA satellites during eclipse seasons. The satellite shapes, area, and optical properties were published for the BLOCK I, II, IIA, IIR [10,11], and IIF satellites [12]. Despite the analytical models having a clear physical interpretation, they can still cause large model errors due to the possible inaccurate satellite structure or optical property. There are two feasible methods to avoid the deviations of these model from the reality. The first one is to adjust optical coefficients with satellite tracking data [13]; for instance, the adjustable boxing-wing modeling considering potential additional acceleration was established using ten years of GPS observations for the GPS BLOCK IIA, IIR, and IIF satellites [8]. The second one is to develop an expression of Fourier series based on the analysis of accelerations caused by SRP, such as the ROCK model series (ROCK-S, ROCK-T), which was first developed by the satellite manufacturer Rockwell International [10,11,14,15], and the GSPM (GPS solar radiation model) developed by JPL (Jet Propulsion Laboratory) [16].
In addition to the models derived from physical-based methods, the empirical SRP models are developed by fitting the long-term orbit estimates regardless of the satellite metadata-for instance, the Empirical CODE Orbit Model (ECOM) developed by left for Orbit Determination in Europe (CODE) [17]. The ECOM was then further refined by reducing the periodic terms in the satellite-Sun direction, that is, the ECOM1 model with five coefficients [18]. The ECOM2 was developed in 2015 to consider the SRP variation caused by different satellite shapes-for instance, the Galileo satellites [19]. Both ECOM1 and ECOM2 are widely used as parameterization models by several ACs [20][21][22]. However, estimating more parameters in ECOM2 can cause an overparameterization and destabilize orbit solutions, and the pure parameterization models may also introduce draconitic errors in GNSS-based geodetic products [23,24]. To minimize the number of SRP parameters and consider the SRP force components that cannot be described by the parameterization model, an a priori SRP model is recommended [25]. Hence, a hybrid strategy of combining the a priori box-wing model and the parameterization model, that is, ECOM1 or ECOM2, has been proved to be beneficial for the orbit quality [8,[26][27][28].
The GPS satellite orbit differences between different ACs have been proved to be largely caused by different approaches of SRP modeling [29]. It was also reported that the performances of different SRP modeling in GPS POD are related to specific satellite types, as most of the BLOCK IIR satellites have different performances [30]. For instance, ECOM1 has better performance than ECOM2 in terms of predicated orbit precision [31]. In addition, an analysis for Chinese BeiDou satellite navigation system (BDS) IGSO and MEO orbits indicated that the box-wing model can remove the βand µ-dependent systematic orbit errors [32]. All the analyses mentioned above have summarized the diversity of GPS orbit precision between using ECOM1 and ECOM2 and the importance of the a priori box-wing model. Despite the above studies, a comprehensive investigation of different SRP models especially when GPS satellites are in eclipse seasons is still lacking.
When a satellite enters the umbra or the penumbra area of the occulting bodies (usually the Earth and the Moon), that is, the eclipse season, a shadow factor is adopted to represent the degrees of Sun occultation. Different ACs adopt different strategies, either in all D direction, which is the satellite-Sun direction, Y direction, which is the solar panel axis direction, or B direction, which completes the right-hand orthogonal sys-tem adopted by GeoForschungsZentrum (GFZ) (Zhiguo Deng, personal communication, 10 February 2020), or only in the D direction adopted by Wuhan University (WHU) (Jin Guo, personal communication, 22 February 2020). Although several studies tried to refine its function model, few studies focus on the relationship between shadow factor and the parameterization model and its impact on eclipsing satellites. The eclipse seasons are quite common for GPS satellites. Figure 1 shows that among the 32 operational GPS satellites in 2019, the average number of eclipsing satellites is around 6 and the maximum number is 12, and a single eclipse event lasts 40 min on average. Moreover, some unmodeled non-SRP forces may exist and will become conspicuous for those satellites during eclipsing season, which cannot be neglected [8,33]. Therefore, it is necessary to investigate the orbit accuracy during eclipse seasons and clarify the proper method of combining shadow factor and the parameterization SRP model. Therefore, the major motivation of this study is to evaluate different combinations of the a priori box-wing model and the parameterization model (that is, ECOM1 or ECOM2), and the application of the shadow factor, focusing on eclipsing satellites. Following this introduction, the SRP models and the shadow function are described in Section 2. After a short description of data processing strategies in Section 3, the results are analyzed in Section 4, where the orbit quality is evaluated by the day boundary discontinuity (DBD), and Earth rotation parameters (ERPs) are evaluated by comparing to the IERS EOP 14C04 product. The main findings are summarized and the future perspectives are given in Section 5.

Methodology
The acceleration of an in-orbit satellite caused by SRP perturbation can be expressed as: where a prior denotes acceleration derived from an a priori model and a para denotes the acceleration represented by a parameterization model to be estimated. In the third IGS reprocessing campaign, the ECOM2, JPL GSPM, or box-wing model combined with ECOM1 is suggested for the GPS satellites except for the BLOCK IIR and IIF, for which the ECOM1 is recommended (http://acc.igs.org/repro3/repro3.html (accessed on 9 May 2021)). For the SRP modeling analysis in this study, we use the box-wing model as the a priori model and ECOM1 or ECOM2 as the parameterization model.
The a priori acceleration a prior can be expressed as [13]: where d denotes the squared ratio between satellite-Sun distance and Earth-Sun distance; A denotes the area of surface; m denotes the mass of the satellite; S 0 denotes the solar flux at 1AU ≈ 1376 W/m 2 ; c denotes the velocity of light in vacuum; α, ρ, and δ denote the fractions of absorbed, specularly reflected, and diffusely scattered photons, respectively; κ denotes the thermal radiation factor (1 for satellite body and 0 for solar panel); e D denotes the unit vector from satellite to Sun; e N denotes the unit normal vector of surface; and θ denotes the angle between e D and e N (cos θ ≥ 0). For the adjustable box-wing model, the optical parameters are updated by tracking data. The acceleration of the parameterization SRP model a para can be expressed as: where µ denotes the satellite's argument of latitude, e D denotes the unit vector in the D direction, e Y denotes the unit vector in the Y direction, and e B denotes the unit vector in the B direction. In ECOM1, the functions D(µ), Y(µ), and B(µ) are represented as a truncated Fourier series [18], where D 0 , Y 0 , and B 0 are constant terms and B c and B s are first-per-revolution terms.
Compared to the ECOM1 model, the twice-per-revolution and fourth-per-revolution terms are added in the D direction in the ECOM2 modeling, while µ is replaced by the angular argument ∆µ = µ − µ s and µ s is the Sun's argument of latitude in satellite orbital plane for straightforward interpretation [19].
Considering the IGS conventional body-fixed frame [34], its relationship with the DYB frame can be presented as: where cos(ε) = cos(∆µ) cos(β), and β denotes the elevation angle of the Sun over the orbital plane. When a satellite enters the eclipse season, a shadow factor f is usually employed to describe the corresponding solar irradiance variation. Adopting the conical shadow model, the shadow factor f can be calculated in the following three cases:

•
Full phase: the Sun is fully visible to the satellite, therefore f = 1; • Penumbra: the satellite is in the penumbra area of the Earth and the Moon and it can receive only part of the solar irradiance from the Sun, therefore [35] (pp. 81-83): where S denotes the total area of the Sun occulted by the Earth and the Moon, and θ Sun the apparent radius of the Sun; • Umbra: the satellite is in the umbra of satellite, therefore f = 0 If the shadow factor is applied in Equation (1) in three directions, which is used by CODE and GFZ, Equation (1) can be rewritten as: If the shadow factor is applied in Equation (1) only in the D direction, which is used by WHU, Equation (1) can be rewritten as: The Y(µ) term usually consists of two parts. The first part is the instrument bias, namely the yaw bias, which is caused by the misalignment of solar panels or solar sensor errors and observed for BLOCK II and IIA satellites [10]. The second one is the thermal radiation, which was found on the Y surfaces (radiated from louvers) of the BLOCK IIR satellites [36]. The thermal radiator is designed for evacuating excessive heat from onboard payloads to space, and as a result, an additional force is generated. However, thermal control information of the GPS satellites is not available to the public, which is hardly to model and cannot be easily included as an a priori SRP model. Ignoring these small forces, for example, using Equation (8) can degrade the orbit accuracy in the eclipse season. An inspection of Equation (9) shows that keeping the Y(µ) active during the Earth's shadow transitions can help to model the thermal radiation during eclipse seasons, along with other possible non-SRP nature forces in Y surfaces. In addition to the Y surfaces, the ±X and −Z surfaces may also undergo thermal radiations. Assuming there is a constant acceleration a −X in the −X surface and it can be represented as a −X sin(ε), 0, and a −X cos(ε) in the D, Y, and B directions, respectively. Similarly, assuming there is a constant acceleration a −Z in the −Z surface and it can be represented as −a −Z cos(ε), 0, and a −Z sin(ε) in the D, Y, and B directions, respectively. Therefore, these constant accelerations corresponding to the B direction can be modeled by the 1 cycle per revolution (cpr) terms in both ECOM1 and ECOM2 because of their 1 cpr terms, whereas the rest part corresponding to the D direction can only be partially modeled in ECOM1 and ECOM2, where ECOM2 has better performance due to its additional 2 and 4 cpr terms.
In addition to the above-mentioned potential thermal radiations, there are other potential nonconservative forces that need to be considered during eclipse seasons-for example, a surface undergoing heating up and cooling down periods during shadow transitions, active thermal control of shadow season, etc. If there are thermal radiation or other unmodeled nonconservative forces in Y or other surface, keeping parameters in the Y and B directions active in the eclipse season can compensate for the forces of non-SRP nature to some extent, and we will continue to prove this inference in Section 4.

Data Processing Strategy
In Section 2, we presented the SRP modeling, the ways of applying the shadow factor, and its relationship with possible unmodeled or inaccurately modeled non-SRP forces. As the ECOM2 has more periodic terms than ECOM1 in the D direction, it may have different behaviors for eclipsing satellites. In this Section, we introduce our analysis methods on evaluating the performance of different SRP models, with an emphasis on satellites in eclipsing seasons.
Combining the strategies of using ECOM1 or ECOM2 as the parameterization model, using no model, the box-wing, and adjustable box-wing model as a priori model, and the application of the shadow factor, i.e., in the D or DYB directions, we have 12 different solutions, as shown in Table 1. For the box-wing model, we use satellite metadata and optical properties provided by the manufacturers [10,11], and for the adjustable box-wing model, we use adjusted metadata and optical properties provided by Duan [8]. From the results listed in Section 4, we infer that each satellite may suffer different magnitudes of unmodeled forces, even for the same block type. The established box-wing model does not contain this term. Considering that the purpose of this work is to evaluate different SRP models and shadow factor, we do not introduce the additional accelerations in either X or Y surfaces provided by Duan et al. [8], which are applied for remaining forces after employing the adjustable box-wing model. Table 1. Twelve cases of GPS daily precise orbit determination solutions using different methods in handling the a priori model, the parameterization model, and the shadow factor. Note that the directions in which shadow factor is applied refers to the parameterization model, and in the a priori model, the shadow factor is always applied. See Equations (8) and (9) for details.

Solution
Empirical SRP Model Shadow Factor A Priori SRP Model The GPS observations from the year of 2017 to 2019 were processed on a daily basis, using the Positioning And Navigation Data Analyst (PANDA) software [37]. More details are given in Table 2. To minimize the possible errors caused by integration, the Adams-Bashforth-Moulton predictor-corrector integration method is switched to the Runge-Kutta-Fehlberg (RKF) integration method during the shadow transition, that is, the period from full phase to umbra or inverse.
Unlike the GLONASS, Galileo, and BeiDou satellites where the SLR tracking observations are available, the GPS BLOCK IIA, IIR, IIF, and III satellites used in this study are not equipped with retroreflector arrays, and thus using SLR as an external orbit validation is not possible [47]. It is also not optimal to use the orbit products form the IGS ACs as a reference for comparison, as their SRP modeling and strategies of applying the shadow factor are different with each other and can affect the orbit product. Therefore, the following quantities are used to assess the different solutions: • Analyses of the estimated SRP parameters: the estimated ECOM parameters should not show systematic biases if the a priori model can describe all important characteristics of the satellite. The correlation between the satellite orbit dynamic parameters can also measure the association of these parameters; • Day boundary discontinuities (DBD) of the satellite orbits: 3D distance between the orbital positions of two consecutive arcs are commonly used to assess the internal consistency of the orbit solution; • Agreement of the Earth rotation parameters (ERPs) to the IERS EOP 14C04 product [45]: the GPS technique provides precise polar motion and LoD estimates, and thus the ERP estimates from different solutions can be compared to the IERS product. We calculated the weighted STD (WSTD) of the differences between the GPS POD solution and the IERS 14C04 product [48] as: where dx i denotes the difference between the ERP parameters and reference values, n denotes the number of records, and σ i and σ re f ,i the formal errors of the estimates and the reference product, respectively.

Results and Discussions
In this Section, we first analyze the orbital parameters in Section 4.1 and the orbit precision of the BLOCK IIR and IIF satellites in Section 4.2, and then present the ERP agreement with the IERS EOP 14C04 product in Section 4.3. For the orbit precision analyses, the BLOCK IIA satellites are not presented because precisely modeling its postshadow recovery can be difficult [49] and only one or two satellites are in operation during the time span. The newly employment BLOCK III satellites are also not included, as the adjustable box-wing model is not available until now and they may undergo outgassing phase at their initial operation period, which will influence the comparison [14].

Estimated ECOM Parameters
To investigate the correlation between ECOM SRP parameters, Figures 2 and 3 show the correlation coefficients in the E2DYB solution among the estimated satellite dynamic parameters (PX, PY, PZ, VX, VY, VZ, D 0 , Y 0 , B 0 , B c , B s , D 2c , D 2s , D 4c , D 4s ) of G067, a BLOCK IIF satellite, and EPRs (XPOLE, YPOLE, XPOLE RATE, YPOLE RATE, LOD) during noneclipsing and eclipsing seasons, respectively. We can see that in both higher and lower β angles, the D 2c term is highly correlated with the D 0 and B c terms, with correlation coefficients larger than 0.9 and 0.6, respectively, and the correlation between D 2s and B s term is close to −1, which means that they are highly correlated. Similar correlation can also be observed for the BLOCK IIR satellites (not shown here). Comparing Figures 2 and 3, we can see that the correlation coefficients are higher in eclipse seasons, i.e., the coefficients between D 2c and PX, B c and VX, D 0 and B 0 at lower β angle are larger than that at the higher β. Consequently, the SRP parameters in the D direction can be hardly separated from the parameters in the Y and B directions, as is the interaction between the different ECOM parameters when satellites are at lower β angle.  Figure 4 presents the time series of estimated Y 0 for all BLOCK IIR and BLOCK IIF satellites in the E1DYB solution. The BLOCK IIR satellites show a larger Y 0 value than the BLOCK IIF satellites, which exhibits as a constant acceleration along the satellite Y-axis, that is, the Y-bias. As the accelerations caused by yaw bias are related to the Sun irradiation, the Y 0 term can be omitted as Equation (8) represented if it is dominated by the yaw bias. The magnitude of perturbations caused by the yaw bias for BLOCK IIR satellites is about 0.15 nm/s 2 [8], which is consistent with the behavior of Galileo Full Operational Capability (FOC) satellites [33,50]. However, most of the estimated Y 0 coefficients for BLOCK IIR satellites are larger than 0.5 nm/s 2 , indicating that there are possible thermal forces in the Y direction. Taking unmodeled non-SRP forces, Equation (9) is more reasonable than Equation (8). It is noteworthy that the highest dispersion is visible for the legacy BLOCK IIR-A and BLOCK IIR-B satellites (the upper two panels in Figure 4), and the largest difference is greater than 0.4 nm/s 2 , whereas the variation of the Y 0 coefficient is less significant for Block IIR-M and BLOCK IIF satellites (the lower two panels in Figure 4). This diversity means the magnitude of thermal forces are different between the BLOCK IIR-A, BLOCK IIR-B, and BLOCK IIR-M satellites.  For a further comparison of Equations (8) and (9), taking SVN G050 (BLOCK IIR) and G073 (BLOCK IIF) as an example, the estimated coefficients in the Y and B directions derived from the E1D and E1DYB solutions are shown in Figure 5. The coefficients of SVN G050 in the B direction are shifted by 6 nm/s 2 for better view. As we can see, the estimated coefficient Y 0 derived in the E1D solution shows less scatter between the phases in and out of eclipse seasons. A similar pattern was also reported when an a priori box-wing model was used (see Figure 15 in [50]). It indicates that BLOCK IIR satellites have a relatively noticeable thermal radiation in Y surfaces. In addition, the Y 0 estimates of BLOCK IIF satellite (SVN G073 in Figure 5) have a visible linear trend during eclipse seasons and are reduced after the β angle switches the sign (see Figure 4 in [8]). This needs further investigation. In contrast to the Y 0 of BLOCK IIR satellites, the trend of B c (for instance, SVN G073) shows larger variation. Furthermore, the standard deviation (STD) values of SRP parameters in the 12 solutions are given in Figure 6. For the ECOM1 solutions (that is, using ECOM1 as the parameterization model) shown in the upper panels of Figure 6, the STD of the B c (marked with black X-cross) is improved when an a priori box-wing model is used or a shadow factor is applied only in the D direction, except for BLOCK IIF satellites, whereas other coefficients have comparable STD values between different solutions. The STD difference of B c between solution E1DYB and solution E1D_ABW is around 0.2 nm/s 2 for BLOCK IIR satellites and BLOCK IIF satellites. As for the ECOM2 solutions (that is, using ECOM2 as parameterization model) shown in the lower panels in Figure 6, the STD values of B c do not show any significant reductions either using box-wing models or adding the shadow factor only in the D direction, that is, between E2DYB and E2D_ABW solution. Compared with the ECOM1 solutions, all SRP coefficients in ECOM2 solutions show lower stability except for the Y o term, and the STD values of the second and fourth order periodic terms are greater than 1 nm/s 2 in the ECOM2 solutions. A noticeable decline of D 4s (blue star) is found when the shadow factor is applied only in the D direction and D 2s shows similar but smaller decline. Summarizing the sensitivity of the GPS BLOCK IIR and IIF satellites to the accelerations in the direction of Y and B, they both suffer from the accelerations in the Y direction, which is independent on β angle, whereas the BLOCK IIF satellites show smaller accelerations. The visible changes during eclipse seasons for the BLOCK IIR satellites may be caused by the active thermal control before and after the middle of an eclipse event and can be removed mostly by considering additional accelerations in the X surfaces. The dependence of the accelerations in the B direction on β angle is caused not only by SRP but also by other nonconservative forces, such as the thermal radiation, as mentioned in Section 2. Applying an a priori box-wing or adjustable box-wing model will improve the stability (STD values) of periodic terms (for instance, B c ) but not for the BLOFK IIF satellites. Compared with using Equation (8), using Equation (9) has a positive effect on the stability of periodic terms, especially for the D 4s term of the BLOCK IIR satellites. Note that nonconservative force modeling is beyond the scope of this study; therefore, we only focus on the performance of SRP modeling combined with shadow factor.

Day Boundary Discontinuity of Satellite Orbits
The mean 3D RMS values of orbit DBD for the BLOCK IIR and IIF satellites in different solutions are shown in Figure 7. The method using ECOM2 as a parameterization model (the E2DYB solution) has slightly improved orbit precision in eclipse seasons compared to that using ECOM1, especially for the BLOCK IIR satellites. Using the box-wing as an a priori model (the E1DYB_BW solution) slightly improves the orbit precision for BLOCK IIR satellites and using the adjustable box-wing as an a priori model (the E1DYB_ABW solution) further improves the orbit quality for both the BLOCK IIR and IIF satellites, in and out of eclipse seasons. By comparing each pair of solutions between the shadow factor applied in the D direction and that in the DYB directions, for instance, between the E1DYB and E1D solutions, and also between the E1DYB_ABW and E1D_ABW solutions, we can see that eclipsing satellites have smaller RMS values of DBD in latter solutions with the shadow factor applied only in the D direction. A slight improvement for noneclipsing satellites can also be observed when the shadow factor is applied in the D direction of ECOM2, which can be explained by the stability improvement of higher periodic terms (D 2s or D 4s ) in Figure 6. For the BLOCK IIR satellites, the orbit precision of E1DYB_ABW and E2DYB_ABW, in which ECOM1 or ECOM2 parameterization model is combined with the a priori adjustable box-wing model and the shadow factor is applied in the D direction, is better than 4 cm for of the noneclipsing satellites, and close to 4 cm for eclipsing satellites. As for the BLOCK IIF satellites, using the adjustable box-wing as an a priori model slightly improves the orbit precision, whereas using the box-wing degrades the orbit precision. Therefore, the solutions using the box-wing model are not discussed in the following analysis.
Taking SVN G050 (BLOCK IIR) and SVN G073 (BLOCK IIF) as an example, Figure 8 presents the radial orbit differences between solutions with the shadow factor applied in the D direction and those in the DYB directions. The along and cross components are not given here because the satellite-Sun direction and radial direction of a satellite are close to collinearity when the satellite enters the eclipse season for lower β angle. As shown in Figure 8, orbit differences between two methods of using the shadow factor are observed at lower β angle, in which the differences based on ECOM2 are more visible. When ∆µ is around 0 • (∆µ + 180 • = ±180 • ), that is, the orbit noon, the orbit differences are negative at lower negative β angle and positive at lower positive β angle. Note that BLOCK IIF satellites show a more visible pattern than BLOCK IIR satellites because of the higher periodic terms in the D direction.
Moreover, the radial orbit differences between using different a priori SRP models are investigated in Figures 9 and 10, which can illustrate the impact of different a priori models on radial orbit. For the GPS SVN050 satellite (BLOCK IIR) shown in Figure 9, using the boxwing or adjustable box-wing as a priori model with the ECOM1 as parameterization model leads to the visible peanutlike distribution of the orbit differences, where the differences can be up to ±6 mm. This pattern is likely owing to the variation of effective illuminated area in the Z direction (pointing to Earth). Large antenna arrays, namely W-sensor, are installed in ±X surfaces of BLOCK IIR satellites, which cannot be modeled by only one constant parameter in the D direction. The differences shown in the left panel of Figure 9 are reduced when using ECOM2 as parameterization model, especially between box-wing model and no a priori model, indicating that the box-wing as a priori model does not contribute that much when ECOM2 is used as parameterization model.
In addition to the BLOCK IIR satellites, the orbit radial differences of the BLOCK IIF satellites (SVN G073) are further investigated in Figure 10. As it is shown, the differences caused by applying a priori box-wing model are largely reduced compared with corresponding results of BLOCK IIR satellites shown in Figure 9. Note that when using the ECOM2 as parameterization model, the visible differences at lower β angle between with and without an a priori adjustable box-wing model may be related to the adjustment of shape value in the adjustable box-wing modeling [8].
Due to fact that current box-wing model has a little improvement on BLOCK IIR satellite orbit precision, even a negative effect on BLOCK IIF satellites, we focus on the evaluation of adjustable box-wing model and the methods of applying the shadow factor. The 3D RMS values of DBD with respect to β angle for eclipsing satellites in different solutions are shown in Figure 11. For BLOCK IIR satellites shown in the upper panel, RMS values of DBD in E1DYB (red dot) and E1DYB_ABW (blue dot) solutions increase with β declining during the eclipse period. Applying the shadow factor only in the D direction removes the rising trend mostly for 4 • < abs(β) < 12 • , that is, the E1D (black dot) and E1D_ABW (green dot) solutions. However, in the E1D solution, the RMS values still increase in the case of lower β angle, that is, abs(β) < 4 • . By using the ECOM2 instead of the ECOM1 and applying the shadow factor in the D direction, the RMS values are reduced to a large extent. Note the relatively larger RMS values of DBD up to 6-8 cm for both BLOCK IIR and BLOCK IIF satellites at the beginning of eclipse events, that is,  ∆µ is the angular argument. Upper: differences between using no a priori SRP model and using the box-wing model; lower: differences between using no a priori SRP model and using the adjustable box-wing model. Left and right panels: using the ECOM1 and ECOM2 as the parameterization model, respectively. For the BLOCK IIF satellites shown in the lower panel of Figure 11, an asymmetrical pattern is observed between negative and positive β angles, which is similar with the SRP coefficients Y 0 variation in Section 4.1 (see Figures 4 and 5). All solutions have small orbit differences with each other when the β angle is positive, including the E1DYB and E1DYB_ABW solutions. The large RMS values are also observed for BLOCK IIF satellites at transition from full phase to umbra or inverse. This phenomenon also needs further investigation.
Further statistics of the RMS values of DBD in the along, cross, and radial components from 2017 to 2019 are summarized in Table 3. The solutions with shadow factor applied in the D direction are generally better that those with shadow factor applied in the DYB directions, especially for the cross and radial components. For the BLOCK IIR satellites in eclipse seasons, the RMS values of DBD are reduced from 2.8 cm and 3.1 cm in the E1DYB solution to 2.1 cm and 2.4 cm in the E1D solution in the cross and radial directions, respectively. For the BLOCK IIF satellites in eclipse seasons, the improvement from the E1DYB to E1D solution is only 0.4 cm in the radial direction, and that in the along and cross directions are neglectable. When an adjustable box-wing model is used, we can see the orbit improvement in and out of eclipse seasons, and the improvement of the BLOCK IIR satellites are larger than that of BLOCK IIF satellites. The orbits are further improved when applying the shadow factor in the D direction and using the adjustable box-wing model as a priori SPR model simultaneously, that is, the E1D_ABW and E2D_ABW solutions. The values of E1D_ABW solution for BLOCK IIR satellites out of eclipse seasons are 2.3 cm, 1.7 cm, and 1.7 cm in the along, cross, and radial directions, respectively, and the corresponding values in eclipse seasons are 2.5 cm, 2.4 cm, and 2.7 cm. Moreover, the DBD of the E1D_ABW solution has a comparable accuracy with the E2D_ABW solution. Table 3. RMS values of the orbit DBD in 2017-2019. We distinguish between statistics in and out of eclipse seasons for both BLOCK IIR and BLOCK IIF satellites. A stands for the along direction, C stands for the cross direction, and R stands for the radial direction. The unit is cm. Note that there is no difference between using no a priori SRP model and using the box-wing model; the solutions using the box-wing model are thus not presented.  Figure 10. Differences of satellite orbits in the radial component between different a priori models for G073 (BLOCK IIF) as a function of β and ∆µ angle in 2019. ∆µ is the angular argument. Upper: differences between using no a priori SRP model and using the box-wing model; lower: differences between using no a priori SRP model and using the adjustable box-wing model. Left and right panels: using ECOM1 and ECOM2 as parameterization model, respectively.

Impact on the Earth Rotation Parameters
In addition to the performance of satellite orbit precision, the ERP components are also sensitive to the satellite orbit dynamic modeling [51]. The ERP component estimations of different solutions are compared to the IERS EOP 14C04 product [48] and the WSTD values are given in Table 4. The different methods of modeling SRP and applying the shadow factor have an insignificant impact on the ERP agreement with the IERS product, and the WSTD of each solution is around 43.6 to 44.7 µas, 23.1 to 24.8 µas, 213.2 to 216.0 µas/day, and 8.7 to 9.8 µs/day for the x-pole, y-pole, x-pole rate, and LoD components, respectively. As for the y-pole rate, however, using the adjustable box-wing as an a priori SRP model clearly improves the precision-for instance, with the WSTD value reduced from 205.3 µas/day in solution E1DYB to 184.9 µas/day in solution E1DYB_ABW. The impact of using ECOM1 or ECOM2 models and that of applying the shadow factor in the D or DYB directions on y-pole rate precision are rather diverse. There is insignificant difference between applying the shadow factor in the D direction and in the DYB directions for the polar motion offset and rate, whereas the LoD always shows an improved precision when the shadow factor applied in the D direction. Moreover, the E1D_ABW solution has the smallest WSTD of the LoD (8.7 µs/day), which is improved by about 8% compared with the E2DYB_ABW solution (9.7 µs/day). Figure 12 further presents the power spectra of the LoD in eight solutions, which demonstrates the impact of the SRP modeling. A clear reduction in the 1.04 cpy and 6.30 cpy signals is observed when applying the shadow factor in the D instead of the DYB directions. The only exception is when using the pure ECOM1 parameterization model without an a priori model, that is, between E1DYB and E1D, where the 1.04 cpy signal is larger in the E1D solution than in the E1DYB solution.

Conclusions
In this study, we investigated the performance of different SRP models, including the box-wing and adjustable box-wing as the a priori model, and the ECOM1 and ECOM2 as the parameterization model in GPS POD solutions. We paid special attention to the orbit behavior during eclipse seasons, focusing on the two methods of handling the shadow factor, that is, applying in either the D or the D, Y, and B directions. We processed the daily GPS POD solution using three years of observations and evaluated the results in terms of the SRP coefficients, the day boundary discontinuity values of orbits, and the ERP agreement to the IERS EOP 14C04.
For the GPS BLOCK IIR satellite, the performance of the SRP coefficient Y 0 between in and out of eclipse seasons is consistent with the shadow factor that is applied in the D direction but rather diverse when applied in the DYB directions. The diversity of the coefficient Y 0 in and out of eclipse seasons is possibly stemming from nonconservative forces in ±Y surfaces, e.g., radiator. A notable discrepancy in the behavior of Y 0 for the BLOCK IIR and BOCK IIF satellites is probably related to the difference in satellite design, thermal control model, etc., which needs further investigation in the future. For the coefficients in the B direction, its dependency on β angle is caused not only by SRP but also by other unmodeled non-SRP forces, which can partly be absorbed by the onceper-revolution terms in ECOM1 and ECOM2. Moreover, applying a priori box-wing or adjustable box-wing model can improve the stability of periodic terms in terms of the standard deviation, especially for the D 2s and D 4s terms in the D direction because of the parameters' correlation, except for the BLOCK IIF satellites using the a priori box-wing model. Among all solutions, the one combining a priori adjustable box-wing with ECOM1 with the shadow factor applied in the D direction (that is, the E1D_AWB) has the best SRP coefficient stability with the smallest STD values.
We also demonstrated that using the a priori adjustable box-wing model and applying the shadow factor in the D instead of the DYB directions both can improve the orbit quality, especially in the cross and radial directions. Compared with the E1DYB solution using no a priori SRP and the shadow factor applied in the DYB directions, the RMS values of orbit DBD in the E1D_ABW solution (using a priori adjustable box-wing model with shadow factor applied in the D direction) are improved by 7.7%, 32.1%, and 35.5% for the BLOCK IIR satellites in noneclipsing seasons in the along, cross, and radial direction, respectively, and the corresponding improvements in the eclipse season are 17.8%, 22.7%, and 26.1%. In the E1D_ABW or E2D_ABW solution, the BLOCK IIR satellites show similar behavior in and out of eclipse seasons. As for the ERP agreement to the IERS EOP 14C04, the LoD is improved when the shadow factor is applied in the D direction instead of the DYB directions. The reason of the deteriorated performance of SRP coefficients and DBD values for the BLOCK IIF satellites when applying box-wing model can be attributed to the fact that the shape and optical properties were collected directly from some documents, which might not be accurate enough.
The analyses in this study provide a verification for the presence of unmodeled thermal effects in GPS satellites [8,33]. Although applying the shadow factor only in the D direction during eclipse seasons can absorb some nonconservative forces, some orbit modeling deficiencies still exist, especially for the BLOCK IIF satellites. A further consideration of unmodeled forces (such as thermal radiator, heating, or cooling of surfaces) in each surface will contribute to establishing a more precise SRP model to improve satellite orbits, especially in eclipse seasons.
the Positioning And Navigation Data Analyst (PANDA) software. The authors also would like to thank the anonymous reviewers for their effort and valuable comments to the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.