Extended forward and inverse modeling of radiation pressure accelerations for LEO satellites

For low Earth orbit (LEO) satellites, activities such as precise orbit determination, gravity field retrieval, and thermospheric density estimation from accelerometry require modeled accelerations due to radiation pressure. To overcome inconsistencies and better understand the propagation of modeling errors into estimates, we here suggest to extend the standard analytical LEO radiation pressure model with emphasis on removing systematic errors in time-dependent radiation data products for the Sun and the Earth. Our extended unified model of Earth radiation pressure accelerations is based on hourly CERES SYN1deg data of the Earth’s outgoing radiation combined with angular distribution models. We apply this approach to the GRACE (Gravity Recovery and Climate Experiment) data. Validations with 1 year of calibrated accelerometer measurements suggest that the proposed model extension reduces RMS fits between 5 and 27%, depending on how measurements were calibrated. In contrast, we find little changes when implementing, e.g., thermal reradiation or anisotropic reflection at the satellite’s surface. The refined model can be adopted to any satellite, but insufficient knowledge of geometry and in particular surface properties remains a limitation. In an inverse approach, we therefore parametrize various combinations of possible systematic errors to investigate estimability and understand correlations of remaining inconsistencies. Using GRACE-A accelerometry data, we solve for corrections of material coefficients and CERES fluxes separately over ocean and land. These results are encouraging and suggest that certain physical radiation pressure model parameters could indeed be determined from satellite accelerometry data.


Introduction
Interactions of photons and gas molecules with the surface of satellites lead to forces, which in turn cause orbit perturbations. Besides the atmospheric drag, the two most important non-gravitational forces acting on a satellite are related to the radiation of the Sun and the Earth. With increasing distance from the Earth's surface, the acceleration due to Earth radiation pressure (ERP) becomes less relevant, whereas the effect of solar radiation pressure (SRP) becomes prevalent.
Accurate modeling of these forces is required for precise orbit determination of GNSS satellites (Fliegel et al. 1992;Rodríguez-Solano et al. 2012;Arnold et al. 2015;Steigenberger et al. 2015;Darugna et al. 2018;Bury et al. 2019), B  Institute of Geodesy and Geoinformation, University of Bonn, Nussallee 17, 53115 Bonn, Germany satellite laser ranging (SLR) satellites (Sośnica et al. 2014;Panzetta et al. 2018), and several Earth observation satellites such as Swarm ), radar altimeter and radar imaging satellites Peter et al. 2017;Hackel et al. 2017). LEO gravity field recovery missions, e.g., Gravity Recovery and Climate Experiment (GRACE; Tapley et al. 2004), CHAllenging Minisatellite Payload (CHAMP; Reigber et al. 2002) or ESA's Gravity field and steady-state Ocean Circulation Explorer (GOCE; Floberghagen et al. 2011), carry space-borne accelerometers to remove non-gravitational forces from measured orbit perturbations. In addition, accurate SRP and ERP force models are required to simulate the total non-gravitational accelerations in case accelerometer measurements are partly missing as at the end of the GRACE mission, or of insufficient quality as with the Swarm mission (e.g., Siemes et al. 2016).
The accelerometer measurements onboard GRACE-D, which is the trailing satellite of the GRACE Follow-On mission launched in May 2018, are of lower quality as expected, and thus, a transplant of accelerometer data from the leading satellite GRACE-C is required (McCullough et al. 2019). Using modeled non-gravitational accelerations at the position of both satellites could lead to an improved accelerometer transplant data (ACT1B) in the future, e.g., when following Behzadpour et al. (2019). Furthermore, the ongoing development of space-borne accelerometers (e.g., Christophe et al. 2018) will likely lead to lower noise levels, and applications, which require separating drag and radiation pressure accelerations, would benefit from refining non-gravitational force models. Modeled ERP and SRP accelerations are also required to derive thermospheric neutral densities and winds from measured accelerometer data (e.g., Sutton et al. 2007;Doornbos et al. 2010;Visser et al. 2019). Another application is the lifetime mission assessment including satellite reentry predictions. Improved force models might also be used in further studies on the Earth's energy imbalance (Hakuba et al. 2018).
The magnitude of SRP in terms of the solar flux is about 1360 W/m 2 at an altitude of 400 km, whereas the Earth's outgoing flux is about 280 W/m 2 , to which the longwave flux contributes about twice as much as the shortwave flux. All fluxes vary with the constellation of Sun, Earth, and satellite including the satellite's altitude above the Earth's surface. The solar flux mainly depends on the solar cycle and the solar rotation with periods of about 11 years and 27 days, respectively. The Earth's outgoing fluxes correspond to the incoming solar flux and thus depend on the solar cycle as well as on the diurnal cycle. The accelerations due to ERP and SRP depend on these fluxes as well as on the satellite's mass, surface area, and plate materials. For GRACE, we found an average acceleration of about 3.7 × 10 −8 m/s 2 for SRP and 1.4 × 10 −8 m/s 2 for ERP during January 2010. In comparison, the average acceleration due to atmospheric drag during the same month (which was close to a solar minimum) was 5.1 × 10 −8 m/s 2 .
Analytical SRP models use only the visible solar flux at the position of the satellite, whereas in reality, the flux and its interaction with the satellite's surface material are frequencydependent. Another simplification is that common shadow functions consider the Earth as a sphere (Montenbruck and Gill 2012). In comparison to SRP, analytical ERP models are based on modeled albedo and emission (Knocke et al. 1988). The interaction of incoming fluxes with the satellite's surface are usually modeled with specular and diffuse reflection as well as absorption using thermo-optical material properties.
However, several aspects in the models are not well covered. For example, in SRP modeling, temporal variations in the solar flux and the solar spectrum are not considered appropriately. Omitting the Earth's flattening and the impact of the atmosphere on the ray of sunlight causes errors in semishadowed regions. Reflections at the satellite's surface are incomplete, since anisotropic reflection is generally ignored and thermal reradiation is not always considered. In common ERP models, the use of albedo and emission data introduces errors that can be avoided nowadays due to the availability of observed fluxes at the top of atmosphere and the consideration of angular dependence of Earth's radiation. Furthermore, thermo-optical material properties and their variation over time are not well known.
For GNSS satellites, analytical SRP models cannot be applied directly, since detailed surface areas and thermooptical material properties are not available. Instead, empirical parameters are commonly estimated Arnold et al. 2015). However, analytical force modeling will become more relevant for Galileo satellites, since optical properties have been published (Bury et al. 2019).
Early analytical ERP and SRP models had been developed to improve the orbit determination of SLR satellites (Rubincam and Weiss 1986). In this context, also the anisotropic reflection at the Earth's surface had been studied (Rubincam et al. 1987). Since the material properties of these spherical satellites represent the largest uncertainty in analytical force modeling for SLR satellites, a radiation pressure coefficient or a scaling factor is commonly estimated within a precise orbit determination (POD;Bloßfeld et al. 2018;Panzetta et al. 2018;Hattori and Otsubo 2019).
In ERP modeling, the Earth's outgoing radiation is usually taken from the Clouds and the Earth's Radiant Energy System (CERES). The CERES instrument consists of a radiometer sensor measuring the radiation at the top of atmosphere, commonly assumed to 20 km, at three spectral channels in slant direction (Wielicki et al. 1996). In the CERES data processing, measured radiation is converted to radiative flux by applying empirical angular distribution models (ADMs; Su et al. 2015a, b). To obtain hourly sampling, the SYN1deg-1hour products used in this paper combine measurements from two CERES instruments with hourly optical data from geostationary satellites (Doelling et al. 2016). In comparison with the monthly available Energy Balanced and Filled (EBAF) data, the SYN1deg products are not constrained to the heat content from ocean reanalysis (Loeb et al. 2018), and thus contain biases, which are expected to introduce systematic errors in ERP modeling. On the other hand, EBAF data rely on a given ocean heat content rate that may not reflect current reality sufficiently.
Our hypothesis is that a careful parametrization of potential systematic errors, based on an extended and consistent forward model for SRP and ERP, and on sensitivity studies of what can be observed with today's space accelerometers, may lead to a new class of empirical radiation pressure models. These models could aid in improving POD, gravity, and thermosphere recovery, but they could also provide clues on systematic errors in radiation data products that we use in SRP and ERP modeling.
In this contribution, we do not address the parametrization of thermosphere density modeling errors. We are aware that contemporary thermosphere models such as NRLMSISE-00 (Picone et al. 2002), JB2008 (Bowman et al. 2008), DTM-2013(Bruinsma 2015, or TIE-GCM (Qian et al. 2014) have been developed based on large data sets but show systematic differences and will likely never be able to predict, e.g., the fast response of the thermosphere to geomagnetic storm events. However, data-assimilating models such as TIE-GCM may achieve the required temporal resolution in the future. On the other hand, it is difficult to model the gas interaction with the satellite's material depending on orbital height. In this paper, we assume that these errors have been sufficiently mitigated such that their effect is random (e.g., fast events) or that errors are simply absorbed in the radiation pressure modeling. Future work will address the problem of joint inversion of thermosphere and radiation pressure errors. We also assume that accelerometers have been calibrated, e.g., within a POD procedure (Van Helleputte et al. 2009;Vielberg et al. 2018).
This study does not address lunar radiation pressure modeling (e.g., Floberghagen et al. 1999). Considering the outgoing radiation of the Moon (Matthews 2008) at an altitude of 400 km from the Earth's surface leads to a flux of about 0.02 W/m 2 . This is less than 0.01% of the Earth's outgoing flux and therefore inconsequential for radiation pressure force modeling for LEO satellites.
In this paper, Sect. 2 presents an extended and consistent forward model for ERP and SRP with a focus on a GRACE-like satellite, i.e., a LEO satellite with flat surfaces. Limitations of conventional standard models are also discussed. In Sect. 3, the results of the forward modeling are presented. In addition, Sect. 4 introduces an inverse model to provide a sensitivity study of potential systematic errors. We summarize the preliminary results from GRACE accelerometery in Sect. 5.

Theory
A satellite's acceleration due to the radiation pressure of the Sun a SRP and the Earth a ERP is in general analytically modeled for a single flat surface by and Equations (1) and (2) result from combining the modified radiation pressure equation in Marshall et al. (1992, p.10) with the wavelength-dependent representation of the flux (e.g., Prölss 2012). Both accelerations depend on the areato-mass ratio A/m. Here, γ denotes the angle between the incident radiation and the normal vector of the surface panel.
In case of ERP, integration over the satellite's field of view Δ at the Earth's surface is required, while considering all surface elements ω. Both the radiation pressure coefficient c R and the radiation pressure P depend on the wavelength λ, and the spectrum of P differs for the Earth ⊕ and the Sun . Additionally, SRP depends on the shadow function ν, which indicates whether the satellite is located in direct sunlight, in shadow or in semi-shadow. ERP and SRP accelerations are commonly computed in an inertial coordinate system and can then be transformed to a satellite body-fixed coordinate system using the satellite's attitude data.
For this study, we consider it prudent to define a 'standard' for ERP and SRP models as a reference for our suggested model extensions. The standard SRP model uses a constant solar flux of 1360.8 W/m 2 (Wild et al. 2013) and considers the Earth's shadow with the assumption of a spherical Earth (Montenbruck and Gill 2012). We refer to the model developed by Knocke et al. (1988) as the standard ERP model, which is based on modeled albedo and emission data. There, the Earth's surface is discretized by 19 elements. In both the ERP and SRP standard models, the interaction of incoming fluxes with the satellite's surface is modeled with specular and diffuse reflection as well as absorption using thermooptical material properties.
In Sects. 2.2 and 2.3, we suggest extensions in radiation pressure forward modeling with a focus on a GRACE-like satellite. We also discuss the approximations commonly made and the consistency of the used data sets. Our final extended model is presented in Sect. 2.4.

Solar radiation pressure
The area A of the satellite's surface panels and the panels' orientation in a satellite frame forms the geometric satellite model. For GRACE, we refer to Bettadpur (2012). While very recent studies Wöske et al. 2019) consider finite element models, here we use the eight-panel model from Bettadpur (2012). Errors in these models are usually not specified, even though they would be helpful for developing a realistic error budget.
The satellite's mass m in Eqs.
(1) and (2) is required to convert the radiation pressure forces to accelerations. During a mission's lifetime, the mass decreases due to the consumption of fuel. Assuming a linear mass loss as in Tapley et al. (2007) can only be seen as a rough approximation. In the case of GRACE, the mass product MAS1B provides mass measurements every 10 s based on (a) thruster firing, and (b) pressure and temperature sensors in the tanks (Bettadpur 2012). The measurement error of (a) is 0.2 kg, whereas the error of (b) is stated 0 for the entire mission in the MAS1B data. In this paper, we apply mass data from tank sensors.
Photon energy depends on the wavelength λ, and the SRP acceleration thus depends on the solar irradiance spectrum. Integrating over all wavelengths of the spectrum yields the total solar irradiance of approximately 1362 W/m 2 (Dewitte and Clerbaux 2017) and 1360.5 W/m 2 during the solar minimum in 2008 (Kopp and Lean 2011), respectively. SPR models commonly use visible wavelengths only (e.g., Sutton et al. 2007;Cerri et al. 2010;Doornbos 2012;Wöske et al. 2019) as input, which is-at least for GNSS satellites-due to the lack of thermal material properties (e.g., Ziebart 2004;Hackel et al. 2017). When the material properties given are in the visible and infrared domain, we can discretize the λintegration in Eq. (1) by an equally weighted summation of visible and infrared light. However, there is usually no information provided on the wavelength band, at which given material properties are strictly valid. We note that a similar integration over all wavelengths of the solar spectrum is also used in thermosphere modeling (e.g., TIE-GCM, Qian et al. 2014) to account for the interactions of thermospheric constituents with the solar irradiance at discrete spectral lines. For this, typically a solar extreme ultraviolet proxy model, e.g., EUVAC (Richards et al. 1994), is applied. Providing the material properties for multiple wavelengths such as ultraviolet radiation would thus allow to discretize the integration over all wavelengths in more detail for future missions and consistent with thermosphere density modeling.
The radiation pressure of the Sun P (λ) is commonly computed from the solar radiation pressure at one astronomical unit (AU), which is related to the radiation pressure at the current position of the satellite using the inverse square law where r ,sat is the distance from the Sun to the satellite, and the position of the Sun can be obtained, e.g., from JPL DE421 ephemerides (Folkner et al. 2008). Solar radiation pressure at 1 AU is commonly approximated by the ratio of the solar constant and the speed of light. But the solar constant is a temporal average of the Sun's energy flux integrated over all wavelengths λ, and since this energy flux varies with the solar rotation and the solar cycle, we replace the solar constant of 1362 W/m 2 (e.g., Dewitte and Clerbaux 2017) with a times series of the daily total solar irradiance (TSI), 1 which varied between 1357.0 and 1362.7 W/m 2 since the year 2000.
This time series has been derived from multiple TSI measurements, e.g., from the Total Irradiance Monitor (TIM) onboard the SOlar Radiation and Climate Experiment (SOURCE). To account for the temporal variation in the Sun's radiation pressure, Eq. (3) then becomes time-dependent The radiation pressure coefficient c R (λ) models the absorption and reflection of incoming photons (here from the Sun ) at the satellite's surface in dependency of the wavelength λ. Commonly, diffuse and specular reflection are considered together with absorption (e.g., Doornbos 2012; Montenbruck et al. 2015) as where s is the normalized direction of incoming radiation and n is the surface normal; the radiation pressure coefficient is zero, if the incident angle γ between s and n is smaller than 90 degrees. In Eq. (5), the thermo-optical properties of the satellite's materials are essential. While in theory for each surface material, the amount of photons that are reflected specularly c s , diffusely c d , and absorbed c a is required for each wavelength, the thermo-optical properties are either assumed to be constant for all wavelengths or separated for visible and infrared radiation, which simplifies the integration over the wavelengths in Eqs. (1) and (2). The material properties of GRACE are provided for visible and infrared radiation within the macro-model (Bettadpur 2012). It is not common to consider changes of the thermo-optical properties, which, however, occur due to (1) the thickness of the coating material (e.g., Silverman 1995 , Table 4.7), (2) the time of UV radiation exposure measured in equivalent solar hours, and (3) the atomic oxygen erosion (e.g., Silverman 1995, p.2.23).
Additionally, some models account for thermal reradiation of the satellite by considering a simple static temperature model (Fliegel et al. 1992;Montenbruck et al. 2015) or by using an advanced model accounting for transient heating and heat conduction (Wöske et al. 2019). Anisotropic or delayed thermal reradiation including the possibly rerouted heat generated by the satellite's electrical components (Rievers et al. 2010) may be relevant, but cannot be considered here. In this paper, we assume simple instantaneous reradiation of heat as in Montenbruck et al. (2015). Rearranging the radiation pressure coefficient [Eq. (5)] and including instantaneous thermal reradiation yield The dependency of the scattered reflection at the satellite's surface on the angle of illumination is rarely considered. Bidirectional reflectance distribution functions (BRDF) account for anisotropy, energy conservation, diffuse, and specular bidirectional reflectance (e.g., Ashikhmin and Shirley 2000). Thus, including a BRDF should make the radiation pressure force models more realistic provided the underlying assumptions are correct. In this study, we apply the Ashikhmin-Shirley BRDF (Ashikhmin and Shirley 2000) as suggested by Wetterer et al. (2014) for radiation pressure force modeling, which is one possible implementation of a BRDF. The two exponential factors in the model (Wetterer et al. 2014) that account for anisotropic reflection are set to 1 because this choice minimizes the difference to the current model [Eq. (6)]. Including the BRDF in terms of correction factors Δ d , Δ s1 , Δ s2 for the thermo-optical material properties leads to an extended radiation pressure coefficient In Eq. (1), the shadow function ν indicates whether the satellite is located in direct sunlight, in shadow, or in semishadow. The most basic shadow function is a cylindrical model that only distinguishes between sunlight and shadow (e.g., Hubaux et al. 2012). However, the penumbra regions cannot be ignored. For example, during January 2010 when GRACE's beta angle was about 37 • the penumbra transition lasted nearly 2 min. Geometrically more advanced are conical models that enable to distinguish between sunlight, shadow, and semi-shadow. These conical models often assume a spherical Earth (Montenbruck and Gill 2012), and sometimes an oblate Earth (Adhya et al. 2004;Srivastava et al. 2014), which is obviously physically more correct. Nevertheless, all geometrical models lack the consideration of the light's absorption, scattering, and refraction in the atmosphere. These atmospheric effects are accounted for in the physical model by Robertson (2015) together with the assumption of an oblate Earth. Due to the complexity and computational effort of the original model, here we apply the curve-fitted model SOLAARS-CF (Robertson 2015). When using the physical shadow function, the penumbra transition of GRACE in July 2010 takes one minute and 20 s, which is more than two times longer than in case of applying the conical model. This causes a difference of 2 × 10 −9 m/s 2 in the modeled acceleration during one penumbra transition.
Additionally, we account for lunar eclipses with the conical model by Montenbruck and Gill (2012), with the assumption of a spherical Moon. Since overlapping of the Earth's and the Moon's shadow is very rare (Zhang et al. 2019), we do not account for these events in this study. We apply the JPL DE421 ephemerides (Folkner et al. 2008) to model the shadow of the Earth and the Moon.
Another shadowing effect is self-shadowing, which is often ignored in satellite force modeling. Mazarico et al. (2009) found that the effect of this phenomenon strongly depends on the shape of the satellite. Löcher and Kusche (2018) apply self-shadowing within an SRP model for the Lunar Reconnaissance Orbiter, and Kenneally and Schaub (2017) apply ray-tracing models to consider self-shadowing, as well as self-reflection, for the Mars Reconnaissance Orbiter. But both orbiters had a complex shape as compared to GRACE, and SRP plays a much larger role in the lunar and Mars force models. Self-shadowing of the GRACE satellites has been considered in Wöske et al. (2019) but was found to cause a very small effect. Therefore, we ignore the selfshadowing of both solar and Earth radiation in this study; however, the influence of self-shadowing and multiple reflections has to be reviewed when transferring the model to other spacecraft as considered, e.g., in List et al. (2015).

Earth radiation pressure
To be consistent with SRP modeling, the satellite's mass, geometry, and the radiation pressure coefficient are considered in ERP modeling here in the same way. In comparison with SRP, ERP accelerations depend on the Earth's outgoing radiation in the satellite's field of view (FOV). Therefore, integration over (a) the satellite's field of view Δ as well as an integration over (b) the wavelengths λ of the Earth's outgoing radiation is necessary [see Eq. (2)].
The integration over the satellite's FOV is commonly discretized by accumulating over each surface element with an individual area. Knocke et al. (1988) suggested a ring-like discretization of the FOV that results in 19 surface elements and is still a standard (Montenbruck and Gill 2012). Other recent publications assume a grid of 2.5 • in longitude and latitude (Rodríguez-Solano et al. 2012). In this paper, the FOV is discretized with a grid of 1 • in longitude and latitude as, e.g., in Visser et al. (2019); Wöske et al. (2019). Currently, a grid of 1 • appears sufficient, since this corresponds to the spatial resolution of the Earth's outgoing radiation data. We compute the satellite's FOV and the area of each surface element of the Earth as in Doornbos (2012).
For each single surface element, integration of the radiation pressure P ⊕ of the Earth over all wavelengths λ radiated from this element is necessary In Eq. (8) (see Eq. (4) in Knocke et al. 1988), α is the angle between the element's surface normal and the vector pointing from the center of this element to the satellite, c is the speed of light, Δω is the discretized Earth surface element dω, i.e., the surface area, and r Sat is the distance between the satellite and the center of the surface element. Additionally, Eq. (8) requires information on the Earth's outgoing shortwave radiance I SW due to the reflectance of the sunlight, as well as the outgoing longwave radiance I LW due to the Earth's thermal emission. The radiance in ERP models is commonly based on albedo and emission (Knocke et al. 1988), where albedo a is modeled as the ratio of outgoing to incoming shortwave fluxes, and emission e is the ratio of outgoing longwave flux to incoming shortwave flux. Then, Eq. (8) becomes (see Eq. (16) in Knocke et al. 1988), where τ is the illumination function assigning whether the surface element is Sun-lit or not and ζ is the incident angle at the Earth's surface element between the surface normal and the incoming solar flux. Finally, F is the solar flux at the position of the Earth, which is assumed to be constant.
Beginning with a second degree zonal harmonic function (Knocke et al. 1988), the representation of albedo and emission has gone through a considerable development due to the availability of the Clouds and the Earth's Radiant Energy System (CERES) measurements at top of atmosphere. However, the data sets used for albedo and emission differ in recent publications.  2019) use hourly SYN1deg data. As mentioned before, unlike the EBAF product, the SYN1deg product is not constrained to the heat content from ocean reanalysis (Loeb et al. 2018) and thus expected to contain a bias of approximately 4 W/m 2 (Johnson et al. 2016). Adding 4 W/m 2 to the global average of the Earth's outgoing longwave or shortwave radiation can change the radial ERP acceleration by about 1% during January 2010. The EBAF product has monthly resolution, whereas the SYN1deg data is available hourly, which means that one no longer requires using an illumination function to consider the shortwave flux only at the dayside.
In the original ERP model, Knocke et al. (1988) use albedo and emission fields due to the lack of outgoing fluxes to model the Earth's diffuse shortwave and longwave radiance [see Eq. (8)]. Nowadays, CERES observations of outgoing fluxes are available and there is no need to use incoming fluxes anymore. We apply the CERES SYN1deg-1Hour Ed4A observed TOA longwave flux F LW and shortwave flux F SW (CERES Science Team 2019). Using these variables, the integrated radiation pressure in Eq. (8) becomes However, the conversion from flux to radiance still misses considering the angular dependence of Earth's radiation. In the CERES data processing, angular distribution models (Su et al. 2015a, b) are considered to provide global fluxes F dependent on the solar zenith angle θ 0 and independent from the instrument's viewing zenith angle θ and relative azimuth angle Φ as See Fig. 1 in Suttles et al. (1988) for the definition of the angles. Equation (11) (see Eq. (2) in Su et al. 2015a) is based on the observed radiance I and the ADM anisotropic factor R depending on the land surface and cloud types. To be consistent, then computing ERP accelerations at the position of a satellite requires the back-projection of CERES fluxes to the angular dependence of Earth's radiation. Then, Eq. (10) in combination with the rearranged Eq. (11) becomes To our knowledge, this is the first study that considers ADMs consistently in ERP modeling. Applying the state-ofthe-art ADMs (Su et al. 2015a, b) would be fully in line with the CERES processing; however, here we use the ADMs developed by Suttles et al. (1988Suttles et al. ( , 1989 because they are much easier to handle. In general, ADMs depend on the viewing zenith angle and on 12 scene types, which consider cloud coverage and land cover type. We use the cloud area fraction provided within the CERES SYN1deg data at the same temporal and spatial resolution as the used fluxes to distinguish between clear sky (0-5% coverage), partly cloudy (5-50% coverage), mostly cloudy (50-95% coverage), and overcast (95-100% coverage). The land cover type used to determine the scene type differentiates between ocean, land, snow, desert, and land-ocean mix. We generate a global land cover map from MODIS/Terra and Aqua Combined Land Cover Type CMG Yearly Global 0.05 Deg V006 (Friedl and Sulla-Menashe 2015). Since our application requires only five land cover types instead of 17, we adopt snow and ocean areas as they stand, whereas barren land with up to 10% vegetation corresponds to desert, and the remaining types correspond to land. Within a majority voting, which is required to obtain the land cover type at a spatial resolution of 1 • (similar to the available fluxes) instead of 0.05 • , we consider 1 • × 1 • grid cells as coastal region, if they contain at least 20% of ocean as well as any other land type (land, desert, snow). The MODIS land cover maps are available yearly since 2001; however, at a spatial resolution of 1 • land cover changes are minor and Fig. 1 Land cover map used to determine the scene type for the angular distribution model in ERP modeling we decided to disregard them. Thus, we use the land cover maps of the year 2010, which is approximately the middle of the lifetime of GRACE. Our final land cover map is presented in Fig. 1. Since the assigned snow regions are permanently covered by snow, information on seasonal snow coverage is taken from the snow/ice percent coverage, which is available within the CERES auxiliary surface data derived from the National Snow and Ice Data Center.
Additionally, ADMs for the longwave flux depend on the season and the colatitude (see Table 2 in Suttles et al. 1989). ADMs for the shortwave flux depend on the relative azimuth angle and the solar zenith angle (see Table 2 in Suttles et al. 1988).

Extended unified model for SRP and ERP
Here, we summarize the extended unified analytical models for SRP and ERP accelerations. Based on Eqs. (1) and (2), we include the modifications of the radiation pressure models as explained in Sects. 2.2 and 2.3.
The extended SRP acceleration acting on a single flat surface of the satellite can be modeled as The SRP acceleration acting on the entire satellite can be obtained by evaluating and summing up Eq. (13) for all surface panels. For GRACE, the extended SRP model differs from the standard model as it includes the dependency on at least two channels of the solar spectrum (visible and infrared) instead of visible wavelengths only. The solar con-stant is replaced by a time series of the total solar irradiance to account for the variability of the solar flux. The reflection at the satellite's surface is modified by accounting for anisotropic reflection (Wetterer et al. 2014) and thermal reradiation . In addition, we apply a shadow function considering processes in the Earth's atmosphere and the Earth as a spheroid (Robertson 2015).
The extended ERP acceleration model reads Total ERP acceleration can be obtained in the same way as for SRP, i.e., by evaluating and summing up Eq. (14) for all surface panels. To be consistent with the SRP model, the reflection at the satellite's surface is implemented as above.
As mentioned before, we modify the radiation data sets in the standard ERP equation by directly using hourly outgoing fluxes at the top of atmosphere instead of albedo and emission maps. We also account for the angular dependence of Earth's radiation using angular distribution models (Suttles et al. 1988(Suttles et al. , 1989. Instead of discretizing the Earth's surface by a small number of elements (Knocke et al. 1988), we use surface elements j with a resolution of 1 • in longitude and latitude.
In this study, the computation of ERP and SRP accelerations is based on the satellite's position taken from precise kinematic orbits (Zehentner and Mayer-Gürr 2016), but any orbit could be used. Accelerations are modeled in an inertial reference frame, where transformation to the satellite reference frame (SRF) is performed using star camera data.

Forward modeling: extended versus standard model
Evaluating the skill of the extended radiation pressure force model with real data appears challenging because this would require perfect accelerometer measurements and a perfect drag model. Therefore, we start with a comparison of the standard radiation pressure force models to different stages of extension as presented in the previous section for GRACE-A. The results in this chapter are represented in the SRF. The modeled SRP and ERP accelerations (Figs. 2, 3) are shown here exemplarily for 3 h during January 1, 2010. During this time period, GRACE went through approximately two orbital revolutions (Fig. 4). In addition to the ground track, Fig. 5 shows the global maps of fluxes, cloud fraction, and snow coverage, which are required in ERP modeling.
In Fig. 2, based on the standard SRP model (red), extensions are added successively until our final SRP model (black). The magnitude in the cross-track and radial directions is up to 4.7 × 10 −8 m/s 2 and 2.5 × 10 −8 m/s 2 in the along-track direction, and the SRP acceleration is zero during eclipse transition. The direction of SRP acceleration in the SRF depends on the orbit of the satellite and its orientation.
Differences between the extended model and the standard model are pronounced when the satellite passes a latitude of ± 51 • , e.g., around 05:18 and 06:05; however, it remains unclear why the BRDF causes these artifacts. We also found that the semi-shadowed regions in the extended model are more than two times longer due to the consideration of the physical shadow function (Robertson 2015).
The ERP acceleration for the same time period (Fig. 3) is largest along the radial direction as expected. The extended ERP acceleration (black) is up to 2.91 × 10 −8 m/s 2 in the radial direction. The different model versions vary significantly, especially when the satellite passes over the Antarctic region, where the extended ERP model captures more details than the other models due to the use of hourly outgoing fluxes. Figure 6 visualizes the variations by zooming into 45 min of the ERP accelerations in the three directions. Large ERP variations in the Antarctic region are clearly related to the strong outgoing shortwave flux due to large albedo.
Subsequently, we show the magnitude of the standard and extended radiation pressure accelerations, i.e., ERP, SRP, ERP+SRP, for GRACE-A during January 2010 with respect to the argument of latitude (Fig. 7). The norm of the extended ERP accelerations with up to 3 × 10 −8 m/s 2 is about two Areas in gray correspond to data gaps times smaller as compared to the standard model; thus, the standard model may result in overestimated accelerations. In addition, the extended ERP model provides more detailed structures due to the increased temporal resolution (hourly) of the radiation data sets. More structures in the SRP model are related to the BRDF. SRP accelerations clearly show the variations in the Earth's shadow in time. At the night side, the total radiation pressure acceleration acting on the satellite is affected by the Earth's longwave radiation only, which happens in the descending orbit during this month. ERP and SRP totally amount to 1 × 10 −8 m/s 2 .
To point out the differences between the model extensions, we compute mean differences and the root mean square differences (RMSD) of modeled accelerations for January 2010 (Fig. 8). The model versions are abbreviated using a combination of five digits that are explained for SRP in Table 1 and for ERP in Table 2. For example, 00000 represents the standard model, whereas 13111 is the extended ERP model and 11111 is the extended SRP model. In summary, we find that introducing ERP and SRP model extensions has a large impact on the radial acceleration, whereas the impact on the along-track accelerations is minor. For ERP, large changes in the radial direction are expected because the radial acceleration is the largest, and for SRP, this depends on the orbital plane orientation. The ERP mean differences and RMSD are largest (1 × 10 −8 m/s 2 in the radial direction) when changing the discretization of the satellite's FOV from 19 surface elements as suggested by Knocke et al. (1988) to a grid of 1 • in latitude and longitude. Considering the outgoing radiation only on a coarse grid can easily lead to over-or underestimated radiation pressure accelerations, since the acceleration strongly depends on the radiation at the position of these surface elements.  Table 1 for SRP  and Table 2 for ERP In SRP modeling, including the instantaneous thermal reradiation model and the simplified anisotropic reflection at the satellite's surface causes the largest mean differences with up to 7 × 10 −9 m/s 2 . Thus, at this point we conclude that these effects deserve further study. Including the physical shadow function causes the smallest mean differences of 1 × 10 −12 m/s 2 and smallest RMSD of 2 × 10 −12 m/s 2 .
Nevertheless, these changes are important when considering accelerations during penumbra transitions.
In a next step, we proceed to validate our extended forward radiation pressure model using real observations. The accelerometer onboard GRACE provides measurements of the total non-gravitational accelerations acting on the satellite. The relation between the observed acceleration a obs and the modeled radiation pressure acceleration is where a drag denotes the acceleration due to atmospheric drag. We are aware that both the drag model (including gassurface interaction) and the calibration of the accelerometer play a pivotal role when comparing measured to modeled accelerations. The choices that we make here may very well affect the outcome of the following experiments. In our view, there are thus two ways to consider this dependency. (1) One applies a single state-of-the-art calibration (without making use of non-gravitational force models) and a reference drag model, which are kept fixed during the following experiments and allow for the comparison of differently modeled radiation pressure accelerations. Then, in an inverse modeling approach, we can study remaining inconsistencies in the radiation pressure model, keeping in mind that they will be affected by any residual miscalibration or calibration instability and/or drag model errors.
(2) On the other hand, we could attempt to combine the accelerometer calibration with the inverse modeling of the radiation pressure accelerations in a joint approach, either implemented one step or iteratively. This, albeit certainly innovative, approach would require very extensive numerical testing and complicated interpretation of estimates. Thus, we stick to option (1), while option (2) remains for future research.
In the following, we model the drag acceleration [in Eq. (15)] as described in Vielberg et al. (2018). The model is the same for all scenarios and accounts for drag and lift forces as in Doornbos (2012) with an energy accommodation coefficient of 0.93. Thermospheric neutral densities are obtained from NRLMSISE-00 (Picone et al. 2002). Neutral winds are not considered.
To avoid the impact of spiking accelerations due to thruster firings, we eliminate accelerations 30 s before and after each thrust based on the GRACE Level 1B thruster activation data. As a consequence, about 25% of the data are excluded from the subsequent analyses.
Due to their measurement principle, accelerometers have to be calibrated before use in scientific applications. We confront modeled accelerations, which were calibrated following two different strategies: Fig. 9 (top) compares modeled accelerations against measurements, which are calibrated using the recommended a priori parameters by Bettadpur (2009)  a quadratic fit of daily bias estimates derived from a precise orbit determination using data until March 2009. The absolute mean differences are at the order of 2 × 10 −8 m/s 2 in the along-track and 3 × 10 −7 m/s 2 in the cross-track directions.

consisting of constant scale factors and biases from
In the radial direction, where the impact of the force model extensions is largest, the mean differences vary around 3 × 10 −9 m/s 2 . Additionally, we compare modeled accelerations to a second set of calibrated accelerometer measurements (Fig. 9, bottom). Here, the calibration parameters (constant scale factor and daily biases) have been obtained within a precise orbit determination procedure following Vielberg et al. (2018, Sect. 3.2) without making use of non-gravitational force models using data between August 2002 and July 2016. In this case, the mean differences between modeled and calibrated accelerations in the along-track direction are at the same order of magnitude as in comparison with the a priori calibration parameters. We here speculate that the mean difference in the along-track direction is caused by a mean bias in the modeled drag acceleration. In the cross-track direction, the mean differences between modeled and calibrated accelerations according to Vielberg et al. (2018) are about 4 × 10 −10 m/s 2 , which is about three orders of magnitude smaller than in the comparison above. Interestingly, the mean difference between modeled and our calibrated accelerations increases by about one order of magnitude when adding the instantaneous thermal reradiation. One could conclude that this thermal reradiation model probably does not represent reality well; however, the increased mean difference might also be related to a less reliable accelerometer calibration in the cross-track direction, since the performance of the accelerometer is less sensitive in this direction (Flury et al. 2008). In addition, we speculate that neglecting horizontal neutral winds in the drag model might explain the discrepancies in the cross-track direction. When comparing modeled accelerations to a priori calibrated accelerometer measurements, the RMS of the differences (RMSD) is largest in the cross-track direction. This is not the case when using our calibrated accelerometer measurements; thus, the a priori calibration parameters in the cross-track direction might be less reliable than in the other directions, which is in line with the larger residuals of the fitted cross-track biases in Bettadpur (2009). The RMSD is largest with about 4 × 10 −8 m/s 2 in the radial direction when comparing modeled accelerations to our calibrated accelerometer data, which might be related to a remaining bias in the calibration procedure. We suppose again that different model extensions do not change the magnitude of the RMS of the differences in the along-track direction due to a possible mean bias in the drag model.
The RMSD reduction emphasizes the impact of each model extension (Fig. 9). The RMSD reduction is here computed for different model extensions RMSD 2 with respect to the standard model RMSD 1 , i.e., the positive results mean that accelerations from the extended model version are closer to calibrated accelerations than the standard model. In the comparison against calibrated accelerations from Vielberg et al. (2018), we found that extended ERP and SRP accelerations have no impact on the modeled along-track acceleration because here the drag dominates. Interestingly, visible and infrared fluxes in the SRP model bring the modeled cross-track accelerations about 5.5% closer to calibrated observations from Vielberg et al. (2018). Including the variation in the solar flux and the physical shadow function has no significant impact on the total acceleration. On the contrary, considering instantaneous thermal reradiation seems to fit more than 39% worse to the calibrated cross-track accelerations than the standard model. Also, we find that including anisotropic reflection at the satellite's surface impairs the fit of the standard model to calibrated accelerations to 78%. Increasing discrepancies when including thermal reradiation and the BRDF does not necessarily mean that these extensions are unsuitable. However, we recall that these discrepancies can be related to neglecting horizontal neutral winds in the drag model, to errors in the accelerometer calibration, or to our very limited knowledge on the actual reflectance properties. Obviously, more research is needed to find develop, calibrate, and test BRDFs for satellite force modeling based on manufacturer information of the surface materials. When using calibrated accelerations following Vielberg et al. (2018), the RMSD reduction in the radial direction is up to 5%, which corresponds to 4 × 10 −9 m/s 2 . In case of using a priori calibrated accelerometer measurements, the radial RMSD reduction is larger with up to 27%, which is equal to 2 × 10 −9 m/s 2 . In fact, this is an expected result, since applying daily biases (Vielberg et al. 2018) brings the measured accelerations closer to reality (and closer to modeled accelerations) than the calibration with fitted apriori biases (Bettadpur 2009). After introducing thermal reradiation and the BRDF, the improvement in the radial accelerations decreases to 2% (and 12%) with respect to the standard model when comparing to calibrated (and recommended a priori calibrated) accelerations. It also appears possible that these large RMSD reductions are related to errors in the radial accelerometer calibration parameters. In other words and as was suggested above already, combining the accelerometer calibration and the inverse radiation pressure modeling in a joint approach would be ideal to overcome remaining discrepancies. (c nm cos (mλ) +s nm sin (mλ)) P nm (cos θ)) 2 Without c dvis zenith c svis i + cs vis 1 i + cs vis 2 i Δt Without c sir , c dir c svis i + cs vis 1 i + cs vis 3 i 1 − cos 2 γ Without c sir , c dir , c dvis zenith c svis i + cs vis 3 i 1 − cos 2 γ FSW3 + FSW4 θ 2 F SW 5 Without c sir , c dir , c dvis --

6
Without c sir , c dir , c dvis , TSI --"Row digit" denotes which errors have been removed from the full set of errors. "First column digit" defines the parametrization of the error in the material coefficients, which is similar for each surface panel i. "Second column digit" defines the parametrization of the error in the CERES shortwave fluxes. Errors in longwave fluxes are parametrized in the same way. Δt denotes a time difference, γ is the incident angle of the flux and the surface normal. c nm and s nm are spherical harmonic coefficients, and O denotes the ocean function depending on longitude λ and latitude θ

Inverse modeling
Our comparisons in the previous section have shown that extending the standard radiation pressure model leads to significant and in any case non-negligible differences. However, predicting accelerations with the extended model, although it formally includes more physical processes, did not always fit the measured accelerations better than when following the standard approach. We suspect one reason for this is that we do not know, and likely will never know beforehand, the physical parameters in the extended model sufficiently well. In addition, the observed accelerations had to be calibrated (see Sect. 3) and we can safely assume that the estimated calibration parameters are neither error-free nor uncorrelated with respect to each other, which then introduces further errors in the calibrated observations. It is tempting to ask whether the extended radiation pressure model can be empirically improved by estimating some of its defining satelliteor radiation source-related parameters.
In this situation, inverse modeling and sensitivity analysis represent a systematic way to analyze the estimability of model parameters and their likely error correlation. This is also expected to improve the understanding of other nuisance parameter estimations, e.g., the estimation of additional parameters in the gravity field recovery that may absorb radiation pressure modeling errors.

Error parametrizations and sensitivity study
In a first step, we defined a larger number of different error parametrizations for the extended radiation pressure model, i.e., we assumed that virtually all physical parameters such as offsets in the material coefficients are unknown. Then, we attempt at inverse estimates for subsets in Table 3 based on the extended analytical SRP and ERP models assuming the (calibrated) GRACE-A accelerometer data records as input during the entire year 2010. Since we found it impossible from one year of data to estimate all parameters simultaneously, as expected, we gradually limit the number of parameters (see the second column in Table 3). At this step, we judged all parametrizations by condition number, rank deficiency, and parameter correlation.
Visualization of the rank deficiency (Fig. 10) shows that only few parametrizations appear estimable in this experiment. We find that the parametrization of the radiation source-related parameters has a larger effect on the estimability than the surface material-related parameters. The solutions with full rank include multiplicative corrections of the radiation data. However, we conclude that both errors related to the Earth's radiation and the satellite's material properties may be included in the inverse modeling.
For further analysis, here we consider only the solution with the smallest condition number. Further studies will have to address the estimability of solutions with different parametrizations, but this would be out of the scope of the present paper. In this final (best-conditioned) parametrization, we estimate scaling factors for the radiation pressure accelerations in the three directions, biases for the specular visible material coefficients, and multiplicative corrections for longwave and shortwave flux over land and over ocean. Figure 11 shows that the spectrum of normalized eigenvalues of the final parametrization drops off as including more parameters. The condition number of the selected solution is 5.46×10 −6 . A jump in the spectrum indicates that the stability of this solution decreases when including this parameter, e.g., the error of the visible material for the inner satellite port panel, which is due to correlations with the error of the inner starboard panel.  Table 3 to interpret the labels. Case 613, denoted with A, is the solution with the best condition number

Covariance analysis and inverse estimates
It is important to note that the correlation matrix of the chosen parametrization (Fig. 12) is based on the extended analytical SRP and ERP force models and GRACE-A accelerometer data during the year 2010 only. During this year, the satellite went through orbital conditions with a beta angle between ±83 • and altitude variations between 449 km and 504 km, with a corresponding variability, e.g., in the satellite's field of view on the Earth's surface. As expected, we find that the inversion would solve errors of the visible material coefficient for front and rear panels as being 94% correlated, i.e., the errors can hardly be separated (Fig. 12). Both coefficients are 96% correlated with the scale in the along-track direction, which is related to the alignment of front and rear panels with the along-track direction. Interestingly, the coefficient of visible material of the nadir pointing panel is 57% correlated with the error in the shortwave flux. The shortwave flux corrections over ocean and over land are correlated by 57%, whereas there is no such correlation for the longwave flux. We conclude from high correlations of remaining systematics in the radiation pressure models that the choice and interpretation of estimated nuisance parameters, e.g., in gravity field recovery, should be reassessed. We recall that high correlations between two parameters indicate that errors in forward force modeling of the one parameter could be easily absorbed in estimated parameters that relate to the other one.
We estimate the radiation pressure model parameters as described in Sect. 4.1 within a least squares adjustment using real data. Again, the used observations are the calibrated accelerometer data of GRACE-A in 2010, where the calibration parameters were obtained within a precise orbit determination procedure (Vielberg et al. 2018). Accelerations 30 s before and after each thrust are not considered. We assume equal weighting in the least squares adjustment, while we recognize that finding a suitable weighting is difficult, since weighting the data based on the accuracy of the accelerometer, which is one order of magnitude less sensitive in the cross-track direction, leads to an extreme down-weighting of the cross-track accelerations such that the system of equations becomes instable.
The estimated parameters (Table 4) must be viewed with some caution; this is to our knowledge the first attempt at improving radiation pressure modeling from GRACE accelerometery data, and further experiments with the optimal data editing, preprocessing, parametrization, weighting and possibly regularization will be required. Nevertheless, they enable a first discussion of remaining systematics in the input data of radiation pressure models. We find that in our experiment the estimated correction of the shortwave flux suggests that the bias in the CERES shortwave fluxes is larger over ocean than over land, which is possibly related to the expected bias in the SYN1deg product since it is not constrained to heat storage from ocean reanalysis (Johnson et al. 2016). However, this cannot be confirmed for the longwave fluxes. The bias corrections in the visible material coefficient  show large variations. The bias corrections are small for the side and zenith panels, which are larger areas coated with solar arrays. Bias corrections larger than one are not physical and are likely related to high correlations in the assumed parametrization, we thus assume they absorb other unmodeled errors, and we speculate that they may also be indicative when estimated jointly within gravity field retrieval. The formal errors of the estimated parameters are also likely too optimistic because the observations are temporally correlated, which is not yet considered in the adjustment.
In an additional experiment, we estimate only either Earth radiation-related parameters or satellite-related (i.e., thermooptical) parameters. We find that an inverse solution with only Earth radiation-related parameters appears stable with maximum parameter correlation of 64% and the estimated multiplicative corrections change only slightly. When we estimate satellite-related parameters only, we find again large correlations and unphysical estimates. It is clear that we cannot estimate material coefficients for all spacecraft plates. We find it possible to stabilize this estimate numerically by fixing the material coefficients of a single plate (i.e., setting corrections to zero), but the solution for the other plates still oscillates. At this point, we conclude that an improved strategy for parametrization of all parameter must still be developed. This should include a suitable weighting and a proper solution for the collinearity problem of estimating the material coefficients, e.g., via a generalized inverse solution that applies an extra minimum condition. In addition, further efforts on improving the error adjustment of radiation pressure models in a joint approach are beneficial for future applications.

Conclusion
We extended radiation pressure models to overcome inconsistencies in analytical standard models. Our main modifications in modeling Earth radiation pressure accelerations are based on hourly CERES SYN1deg fields of the Earth's outgoing radiation combined with angular distribution models. Validations with one year of calibrated GRACE accelerometer measurements according to Bettadpur (2009) andVielberg et al. (2018) suggest that the proposed model extensions (without thermal reradiation and BRDFs) decrease RMS fits to 5% and 27%, respectively. We conclude that the extended forward model for ERP and SRP accelerations has reduced systematic errors in the analytical models.
Our validation using GRACE accelerometer measurements revealed that including physically more detailed variations in the solar flux and directly using the Earth's outgoing fluxes improves the radiation pressure force model. However, other modifications such as including instantaneous thermal reradiation and anisotropic reflection at the satellite's surface did not lead to improvements.
Instead of further improving the extended radiation pressure model by trial and error, the inverse modeling is a systematic procedure to study remaining inconsistencies. It turns out that the estimability depends more on the parametrization of the radiation source-related parameters than on the material-related parameters. For one year of GRACE radiation pressure accelerations, we found that estimates for errors in the Earth's outgoing shortwave radiation are correlated up to 57% with the visible material coefficient of the satellite's nadir panel. From high correlations of remaining systematic errors, we conclude that the choice and interpretation of estimated nuisance parameters, e.g., in gravity field recovery, should be reassessed.
Conclusions drawn from modeling radiation pressure accelerations for GRACE can also be adopted to other missions as well. A prerequisite is, however, the availability of the satellite's geometry and material surface properties. Applying our extended model on a finite element macromodel is also possible and expected to improve radiation pressure accelerations for GRACE as well.
Further improvement in ERP and SRP accelerations requires access to accurate satellite geometry. In addition, material surface properties for the whole spectrum would allow the consideration of the incoming fluxes at different wavelength bands, which is expected to bring accelerometer observations and models closer together. Long-duration experiments will be beneficial to quantify changes of the material properties related to the solar exposure time and the interaction with atmospheric particles. The extended analytical radiation pressure model might be used to derive empirical model parameters for satellites without proper geometry and material information in the future.