Bandwidth correction of Swarm GPS carrier phase observations for improved orbit and gravity field determination

Gravity fields derived from GPS tracking of the three Swarm satellites have shown artifacts near the geomagnetic equator, where the carrier phase tracking on the L2 frequency is unable to follow rapid ionospheric path delay changes due to a limited tracking loop bandwidth of only 0.25 Hz in the early years of the mission. Based on the knowledge of the loop filter design, an analytical approach is developed to recover the original L2 signal from the observed carrier phase through inversion of the loop transfer function. Precise orbit determination and gravity field solutions are used to assess the quality of the correction. We show that the a posteriori RMS of the ionosphere-free GPS phase observations for a reduced-dynamic orbit determination can be reduced from 3 to 2 mm while keeping up to 7% more data in the outlier screening compared to uncorrected observations. We also show that artifacts in the kinematic orbit and gravity field solution near the geomagnetic equator can be substantially reduced. The analytical correction is able to mitigate the equatorial artifacts. However, the analytical correction is not as successful compared to the down-weighting of problematic GPS data used in earlier studies. In contrast to the weighting approaches, up to 9–10% more kinematic positions can be retained for the heavily disturbed month March 2015 and also stronger signals for gravity field estimation in the equatorial regions are obtained, as can be seen in the reduced error degree variances of the gravity field estimation. The presented approach may also be applied to other low earth orbit missions, provided that the GPS receivers offer a sufficiently high data rate compared to the tracking loop bandwidth, and provided that the basic loop-filter parameters are known.


Introduction
ESA's three satellite mission Swarm was launched in November 2013 (Friis-Christensen et al. 2008). The satellites were placed in polar low Earth orbits with initial altitudes of 480 km (Swarm A, C) and 530 km (Swarm B) after the commissioning phase. The three satellites are equipped with geodetic-grade dual-frequency GPS receivers provided by Rüstungs Unternehmen AG (RUAG).
Due to the gap between the dedicated earth gravity field missions GRACE (Gravity Recovery And Climate Experiment; Tapley et al. 2004) and GRACE-Follow On, Swarm became a gap filler to provide monthly gravity field solutions (Lück et al. 2018). Especially in the first months of the mission, when the solar activity was relatively large, see Fig. 1, and during evening local times, artifacts in the gravity field became visible in Swarm gravity fields around the geomagnetic equator (Jäggi et al. 2016).
The GPS phase measurements leading to these artifacts can be identified using time-derivatives of the L1-L2 carrier phase difference. This difference is directly related to the slant TEC, i.e., the total electron content along the line of sight (Jäggi et al. 2016;Schreiter et al. 2019). Similar artifacts have already been observed for the GOCE (Gravity field and steady-state Ocean Circulation Explorer) Mission (Drinkwater et al. 2007;Jäggi et al. 2015), which used a different type of GPS receiver but identical chipset and tracking technique for semi-codeless tracking of the encrypted P(Y)-code.
After identifying bandwidth limitations of the carrier tracking as a likely cause of systematic measurement errors, the L1 and L2 phase locked loop (PLL) bandwidths were increased several times on the various satellites of the Swarm mission (Table 1; ESA 2015aESA , b, 2016. These updates had a positive impact on both precise orbit determination (POD) and gravity field determination as shown by van den IJssel et al. (2016) and Dahle et al. (2017) but could only benefit observations collected after their implementation. In a first effort, heuristic down-weighting schemes were applied to mitigate the impact of bad carrier phase observations during rapid changes in the L1-L2 carrier phase difference (Schreiter et al. 2019). To cope with this situation, the present study aims to develop a method for correcting carrier phase observations affected by bandwidth limitations. Based on a model of the loop filter provided by the receiver manufacturer, we assess the impact of rapid slant-TEC changes, and propose a method to recover the true L2 signal and to correct the L2 carrier phase observations. Eventually, we investigate the impact of these corrections on orbit and gravity field quality.
We choose March 2015 and August 2015 for our tests to cover the narrowest L2 loop bandwidths (Table 1) with the most significant loop-related tracking errors and also to allow for direct comparison of different loop settings. Since all three Swarm satellites had the same loop settings in March 2015, we focused on Swarm A in this month. This satellite is part of the lower pair and thus affected by a larger slant-TEC than Swarm B. In August 2015, the L2 bandwidth for Swarm C was already updated to 0.5 Hz, but for Swarm A it was still at 0.25 Hz. Since Swarm A and Swarm C fly in close formation, they observe almost the same slant-TEC, which allows for a direct comparison of the tracking performance. For the computation of the orbits and phase residuals as well as the gravity field solutions, the development version of the Bernese GNSS software v5.3 was used (Dach et al. 2015).

Models
This section discusses the origin of bandwidth-related phase tracking errors in the Swarm GPS receiver and provides the mathematical models for their correction. The loop filter response is analyzed in frequency-domain and a loop filter specific transfer function estimated. Eventually the transfer function is used to invert the loop filter.

L2 tracking model and correction
Apart from constant offsets, the L1 and L2 carrier range 1 , 2 may be split into the sum of a frequency-independent geometry term g(t) and a frequency-dependent  ionospheric term, denoted as I(t) for the L1 carrier frequency f 1 . The impact of the ionospheric contribution on 2 may be expressed as f 2 1 ∕f 2 2 I(t) in a first order approximation neglecting the higher-order ionospheric terms: In a time domain representation, the measured phase ̂ (t) is obtained by convolution of the input phase (t) with the loop specific transfer function H(t): Alternatively, the Fourier transform Φ (f ) of the measured carrier phase range is given by the product of the frequency-domain transfer function H(f ) and the Fourier transform of the true carrier range Φ(f ) . In case of L1 carrier phase tracking the tracking loop bandwidth is considered to be sufficiently high for both, the geometric and the ionospheric variation. Accordingly, the measured L1 phase range closely matches the true range: For the much smaller L2 bandwidth, this assumption does not hold. To be insensitive to geometry-related signal dynamics, the L2 PLL is therefore aided using the L1 carrier rate. The L2 carrier phase tracking may thus be described as which expands to The L2 tracking error is therefore given by To recover the true L2 carrier phase measurement, we apply H −1 2 (t) to (6) and obtain with gf = 2 − 1 denoting the geometry-free L1/L2 carrier phase combination. The term in square brackets describes the error of the L2 carrier phase and can be (3)(t) = H(t) * (t).
evaluated based on knowledge of the loop transfer function and the measured L1-L2 carrier phase difference. The inversion is best performed in the frequency domain, where the inverse transfer function is given as H −1 2 (f ) = 1∕H 2 (f ) . It should be mentioned that a GPS carrier phase measurement as provided, for example, in RINEX (Receiver Independent Exchange format; IGS 2019) observation data files is not the direct output of the tracking process. Instead, it is formed from the output of a numerically controlled oscillator (NCO), which is controlled by the actual tracking loop to follow a down-converted version of the received GPS signal. However, due to linearity of the required transformation and the properties of higher-order loop filters, Eq. (9) is not only valid for the measured NCO phase but can likewise be applied to RINEX carrier phase measurements.

Swarm tracking loop
According to personal information obtained from the manufacturer, the Swarm GPS-receivers use a Digital Phase Lock Loop (DPLL) of order 3 with rate-only feedback. An integration interval T = 0.01 s is used for the L1 carrier tracking, which is based on correlation with the open C/A-code signal. In contrast, a value of T = 0.1 s is used for the L2 phase measurement to compensate the incerased noise of the semi-codeless P(Y)code tracking. Also, the PLL bandwidth for L1 is significantly wider than that for L2 (10-15 Hz vs. 0.25-Hz). By aiding the L2 tracking with the phase-rate from the L1 tracking, the capability to follow rapid signal variations associated with the dynamical motion of a satellite in low earth orbit (LEO) can, nevertheless, be retained despite the low L2 bandwidth. The tracking loop implementation involves a one-step delay between estimation and use of the phase rate and is described in Fig. 2.
The formulation itself is given by Stephens and Thomas (1995). We assume an ideal phase extractor. The n-th model residual phase ∼ n is given by where n is the measured phase, and ̂ n is the model phase. In case of a rate-only feedback loop, the (n + 1)-th model phase is given by where T is the integration time, and the dot denotes the first time derivative. With a computation delay of one update interval ̇n +1 T can be obtained using (10) n = n −̂n, The coefficients K 1 , K 2 , and K 3 characterize the tracking loop properties, such as bandwidth, dampening, etc. For the Swarm GPS receiver, the coefficients together with other parameters are given in Table 2.

Continuous-update approximation
In the limiting case of infinitesimally small update intervals T , the discrete-update (DU) loop may be described by a continuous-update (CU) loop with a rational transfer function of order N = 3 in the Fourier domain, where j = K j ∕T j and s = 2i f (Stephens and Thomas 1995). This leads to the CU loop bandwidth Following Ward et al. (2006), this may also be expressed as are bandwidth-independent filter constants that determine the damping and overshoot of the output signal in response to a step change of the input signal.
natural frequency of the filter and determines the filter bandwidth for given a and b. Typical values of these parameters are a = 1.1 and b = 2.4 (Ward et al. 2006), whereas values of 2.5-3 (Table 2) apply for the Swarm PLL in accordance with the choice of supercritical damping. Making use of (13) and the known loop parameters, the signal ĝ(t) at the output of the tracking loop can be computed for a given input signal g(t) by means of the Fouriertransform F using the relation Given H(s) ≠ 0∀s , this may be inverted, such that Use of this relation along with the Fourier-domain representation of the transfer function in (13) provides a computationally convenient and effective way of convoluting a known output signal with the inverse loop transfer function. It avoids the complexity of inverting the discrete-update transfer function in the time domain and will be used in this study for recovering the L2 tracking error based on (9).

Empirical transfer function
While convenient to use, the continuous-update approximation does not provide a realistic model of the actual loop

Fig. 2
Simplified tracking processor with loop-filter, adapted from Thomas (1998) Stephens and Thomas (1995) The values apply for rate-only numerically controlled oscillator (NCO) updates, super-critical damping, and a one-step computational delay. For comparison, continuous-update loop coefficients 0 , a, and b as defined in Ward et al. (2006) are provided for the respective update intervals T behavior if the product of the bandwidth B and the integration time T violates the condition BT ≪ 1 (in Stephens and Thomas (1995): BT ≤ 0.02 ) and if the loop-filter has a computation delay as it is the case for the Swarm GPS receivers. Notable deviations from the true transfer function can, in particular, be noted in the phase shift at frequencies above the characteristic frequency (see Figs. 3 and 4). Here H emp denotes the empirical transfer function derived from an artificial signal, H CU is the continuous update approximation and H fit the approximated transfer function, see (18). To cope with this limitation, we evaluated the transfer function of the discrete-update loop implementation according to Fig. 2 for a signal covering the frequency range of interest. The amplitude and phase response in the frequency domain were then approximated by a rational transfer function with an order 6 in the denominator and order 4 in the numerator (4/6-order), that exceeds the actual loop order, and coefficients a i and b i adjusted to minimize the deviation from the actual transfer function. The specific form and order of the transfer function result from Aguirre and Hurd (1984) for a third of DPLL with computation delay. As shown in Figs. 3 and 4 for the example of the 0.25 Hz bandwidth, a good representation of the DU transfer function for frequencies of 0.001-1 Hz is obtained with a 4/6-order rational transfer function. Independent approximations of this order were determined for all relevant L2 loop settings in Table 2 and are used in the subsequent sections for correction of the corresponding phase measurements.

Tests on synthetic data
To evaluate the tracking loop performance, we use a constant-zero signal with a cosine-shaped pulse of 0.1 Hz frequency and 2 m peak-to-peak amplitude starting at T = 10 s. In addition, white noise with a σ = 1 cm is added. Figure 5 shows the output of the discrete-update tracking loop at T = 0.1s for different bandwidths as used by the Swarm receivers. As expected, Higher bandwidths result in a faster response. However, all loops overestimate the peak and need a few seconds time to follow the zero signal again. In Fig. 5 (bottom), the differences to the original signal are shown. The difference for the 0.25 Hz loop shows the largest amplitude, but the other bandwidths also show clear signatures. This generates systematic biases in the observations provided by the loop when encountering fast variations of the input signal. When applying the inverse of the empirical transfer function derived in the previous section to the output signal with a sampling interval of 1 s, i.e., 10 • T , one may recognize a close match between the original signal and the recovered signal (Fig. 6). Thus, the test demonstrates the effective use of a rational approximation of the transfer function for correcting bandwidth-related tracking errors of a discreteupdate PLL, even if the measurement sampling rate is notably lower than the actual loop update rate but still comparable to the bandwidth and signal frequency.

Application to real data
For the correction of Swarm GPS observations based on (9), we make use of the geometry-free linear combination Fig. 3 Amplitudes of the discrete-update transfer function, a 4/6order approximation, and the continuous-update approximation for the Swarm L2 tracking loop with B = 0.25Hz and T = 0.1s Fig. 4 Phase shifts of the discrete-update transfer function, compared to a 4/6-order approximation, and the continuous approximation for the Swarm L2 tracking loop with B = 0.25Hz and T = 0.1s gf =̂ 2 −̂ 1 , where ̂ 1 and ̂ 2 denote the observed carrier phase range on the L1 and L2 frequencies as obtained from the RINEX observation files at a 1 s sampling. No dedicated cycle slip detection and correction is applied, since most cycle slips are already corrected in the Swarm Level 1B data (NSI 2019). However, the observation arcs were split if gaps of more than 1.5 s or if | � gf (t k ) − � gf (t k+1 )|∕Δt > 1m∕s . For each of the obtained continuous phase arcs obtained this way, we will perform the inversion independently.
As discussed above, the actual discrete-update transfer function for the loop settings applicable on the day of interest was approximated by a rational function H fit (s) in the frequency domain (cf. (18)). The convolution of ̂ gf with the inverse transfer function in (9) was then performed by multiplication of H fit (s) with the Fourier transform of ̂ gf and back-transformation into the time domain. For application on the real phase observations, which are collected at equidistant epochs t k , k = 0, … , N , we use the Fast Fourier Transform (FFT). The phase arcs are not periodic but are assumed to be periodic by the FFT. According to Fig. 5, the main response of the loop Filter takes place within 30 s after a signal. Therefore, discontinuities at the edges need to be avoided. For that purpose, we use a polynomial fit of degree 1 to the first and last 20 s of the phase arc and extrapolate the signal 60 s at the beginning and the end. To obtain a smooth transition between the polynomial and the original signal, we blend in the first and last 10 s using a linear weighting. Eventually, we also detrend the extended signal using a linear function such that no large step occurs in a periodic continuation of the signal. This transformation does not affect the loopfilter response and thus the result of (9), since a tracking loop of order three has no tracking error due to phase velocity (Ward et al. 2006). Note that the extended data vector is only used to avoid edge effects in the convolution with the inverse transfer function. Data outside the original interval are discarded after computing the L2 correction.

Results
The precise orbit determination was performed using the development Version Bernese GNSS Software 5.3. As external products, the final GPS orbits and the high-rate 5 s GPS satellite clock corrections are used, which are provided by the Center Of Orbit Determination Europe (CODE; Dach et al. 2017). The Swarm GPS RINEX files (baseline 0401) and attitude (baseline 0401) are provided by ESA. Further inputs for the reduced dynamic orbit determination include the gravity field model EGM2008 (Pavlis et al. 2012) and the ocean tide model FES 2004 (Lyard et al. 2006). A set  of constrained piecewise constant accelerations in radial, along-track, and cross-track direction is estimated at sixminute batch intervals to account for model deficiencies and non-gravitational forces. The kinematic positions are computed from undifferenced GPS carrier phase observations. A first estimate for the corrections needed is the ionosphere-free (IF) phase residuals to a reduced dynamic (RD) orbit and the associated receiver clock solution. Since the RD orbit shows a higher dynamical stiffness, it better represents the assumed "true" position of the satellite than the kinematic positions. At the locations where artifacts occur, the IF phase residuals may become large (Schreiter et al. 2019), but due to the estimation of epoch-wise receiver clock offsets, they have an approximately epoch-wise zero mean. Figure 7 (top) shows an example of IF phase residuals for the GPS satellite G01 and the associated measurement corrections. For comparing the corrections to the IF phase residuals, the corrections must be scaled by approx. 1.546 due to the pre-factor of 2 in the IF linear combination. At the locations where the IF phase residuals are getting large, the corrections show very similar behavior. However, earlier (around second 12,200 s), there is an apparently opposite correction for G01, which is caused by another GPS satellite (G23). The epoch wise phase residuals are coupled by the epoch-wise receiver clock estimate. As G23 is highly affected, the receiver clock estimate is also affected, which results in a different estimated range for G01.
After applying the corrections, the ionosphere-free phase residuals with respect to a reduced dynamic orbit are notably reduced after applying the corrections. In Fig. 8, one equator-crossing pass is displayed with the black line indicating the equatorial crossing. Around 48,000 s and 50,500 s are the polar regions, which typically show larger residuals. Near second 49,000 s and second 49,500 s larger spikes are visible. These spikes then lead to systematic differences when comparing a reduced dynamic orbit to kinematic positions. It may be recognized that after applying the corrections to the L2 phase observable, the residuals become smaller in the polar regions, but also the spikes disappear to some extent. However, the negative part of the first spike is still present. It is the beginning of a phase arc, where no reliable correction could be performed.

Orbit quality
To evaluate the impact on orbit fitting level, we compare the number of observations used for the final orbit determination, the post-fit root-mean-square (RMS) for both the reduced dynamic and the kinematic orbit and the number of ambiguities set for each scenario, see Figs. 9, 10, 11 and 12. It can be observed that more observations could be used for all days if the corrections to the L2 phase observable are applied. This is mostly due to smaller IF-residuals, which  reduces the number of observations that are rejected in the phase screening. Also, this leads to fewer gaps and, in turn, allows for a lower number of ambiguity parameters (Fig. 12). For March 2015, up to 7% more observations could be retained. For August 2015 still up to 1% more observations are used. Figure 9 (bottom) also shows an increased number of useful observations for Swarm C that uses a twotimes wider bandwidth in August 2015 than Swarm A. This increase relates to the fact that a fixed threshold is used on ionosphere-fee phase residuals in the preprocessing to screen for bad observations in all data sets. As a result of the wider bandwidth of Swarm C, systematic tracking errors and the magnitude of phase residuals are reduced compared to Swarm A. Accordingly, a larger number of observations are accepted and used for POD. The increased number of available observations due to the tracking loop updates were already observed in van den IJssel et al. (2016).
The post fit RMS could be improved from 3 mm RMS down to 2 mm RMS for the reduced dynamic orbit for March 2015. Even the very disturbed days (17.3.2015-19.3.2015) show an improvement, however, the RMS is still highly above typical levels due to a geomagnetic storm taking place on 17.3.2015 with Kp-indices up to 8− (GFZ 2019). Also, for the kinematic positioning, an improvement of the carrier phase residuals from 2.1 mm down to 1.6 mm can be observed. For August 2015 this difference is less pronounced, but still the RMS regarding the reduced dynamic orbits could be reduced from 1.9 mm to 1.5 mm for Swarm A and from 1.6 to 1.4 mm for Swarm C. For the kinematic positioning a reduction from 1.7 mm down to 1.3 mm, Swarm A, and 1.5 mm to 1.25 mm, Swarm C was achieved. Not a single day was degraded.
For an external assessment of the orbit quality, we make use of Satellite Laser Ranging (SLR) validation Arnold et al.  (2018). Even though the orbits were processed using the ITRF2008-compatible reference frame realization IGb08, we use the SRLF2014 instead of the SLRF2008. This approach was proven to be beneficial in Arnold et al. (2018) because of information for post seismic deformation and improved station coordinates. We use a subset of 12 of the available SLR stations: 7090, 7105, 7119, 7501, 7810, 7825, 7827, 7839, 7840, 7841, 7941, and 8834. Among these stations are Herstmonceux, Graz, Greenbelt, Mount Stromlo, Yarragadee, and Zimmerwald, which are known for particularly high quality and amount of SLR observations. An outlier threshold of 200 mm is used. A selection of stations near the geomagnetic equator could not be used due to the very limited number of active SLR stations in that region. For the reduced dynamic orbits in March 2015, the results are given in Tables 3 and 4. Here, both the mean value and standard deviation are reduced by 0.5-1.0 mm, and the RMS SLR residuals decrease by about 1.2 mm. For the kinematic orbit standard deviation and RMS also improves in the corrected scenario, the offset drops by 0.9 mm. The standard deviation and the RMS are reduced by approximately 3.6 mm. For August 2015, in Tables 5 and 6, no large differences can be observed in the reduced dynamic scenario. However, for the kinematic orbits, the mean offset is reduced by 0.2 mm   to 0.3 mm. Only minor improvements at the sub-millimeter scale can be observed in RMS and standard deviation. For August 2015, the ionospheric activity was much lower, see Fig. 1, and the updates of the L2 PLL were already performed (Table 1). Also, the lack of SLR stations near the geomagnetic equator weakens the capability to validate possible orbit improvements in that region.

Gravity field solutions
The gravity field solutions were computed using the celestial mechanics approach (Beutler et al. 2010). Further, we follow the procedure outlined in Jäggi et al. (2016). In the celestial mechanics approach, we first fit an a priori orbit to the kinematic positions considering the EGM2008 gravity field model and the ocean tide model FES 2004. Daily arcs are computed with 15 min. empirical piecewise constant accelerations to compensate for deficiencies in the orbit model. Normal equations for the orbit and gravity field parameters are set up on a daily basis, using again the kinematic positions as pseudo observations. We set up the normal equations for each day and pre-eliminated the orbital parameters and empirical accelerations. Then, we stack the normal equations of one month and solve for the gravity field coefficients. To evaluate the gravity field solutions, we compare to the monthly JPL-RL06 GRACE gravity field solution (Bettadpur 2018;GRACE 2018). First, we perform a visual inspection of the resulting geoid height differences, smoothed using a Gaussian filter with a radius of 400 km, checking if the equatorial artifact is mitigated and if other artifacts occur (see Fig. 13). The correction scenario is capable of reducing the equatorial artifact to a limited extent. We also compare to a solution obtained in a previous study using weighting of observations (Schreiter et al. 2019). Still, the equatorial artifact is least visible in the weighting solution. However, the noise patterns in the polar regions are less pronounced when correcting L2 measurements based on the inverse loop transfer function. In the weighting approach, observations with loop-related tracking errors are notably down-weighted, which leads to large co-variances in the kinematic positions and weakens their impact on the gravity field solution. Using the kinematic positions based on the corrected observations, these positions still have a similar weight in the gravity field solution as in the original one, since the co-variance given for the positions mostly represents the geometry (Jäggi et al. 2011). In the difference and error degree amplitudes, Fig. 14, the corrected observations lead to a reduction of the difference w.r.t. the reference gravity fields and show the smallest error degree variances (dashed line). However, still, the difference amplitudes are significantly smaller using the weighting strategy. Due to insignificant differences, this plot is not displayed for August 2015.
The solution based on L2 phase corrections and the solution based on data weighting show an improvement compared to the original solution (Fig. 14). Formal errors are smaller at all degrees when using the corrected L2 data as compared to the original data due to a larger number of accepted observations. For the actual errors, the benefit of the L2 correction can primarily be seen in higher degrees (> 10). However, this may be expected, as the correction mostly affects frequencies near 1/30 Hz, which corresponds to a few hundred kilometers in spatial resolution and is therefore invisible in the low degrees, which cover larger scales. In the very low degrees, the weighting scenario is slightly better, whereas in the higher degrees, the L2 correction scenario and the weighting solution become comparable. Comparing the RMS of geoid height differences weighted with the cosine of the latitude (wRMS) for March 2015, see Table 7, again the weighting solution shows the lowest value with 12.13 mm compared to the monthly GRACE JPL-RL06 gravity field solution. Still, the best fit with respect to the kinematic positions and the maximum number of kinematic positions used are obtained in the correction scenario. When comparing the results for August 2015, see Table 8, one may see that the wRMS of 8.8 mm, almost matches the value from the weighted solution, 8.5 mm, for Swarm A, but for Swarm C the weighted RMS is 8.8 mm for the weighted solution and 8.6 mm for the unmasked correction. In the  wSTD over the oceans, the corrected scenario shows the smallest value. The number of positions used is similar in all cases, but the best fit with respect to the kinematic positions could be obtained using the correction.

Summary and Conclusions
The limited bandwidth for L2 carrier phase tracking in the Swarm GPS receiver is responsible for increased carrier phase errors occurring near the geomagnetic equator and related artifacts in kinematic position solutions and gravity field models. Whenever the tracked GPS signals are subject to rapidly changing ionospheric path delays, the loop filter introduces systematic biases in the measured L2 phase, whereas the L1 phase is essentially unaffected due to much higher bandwidth. Based on the L2 PLL design knowledge, the loop response can be modeled, and the L2 tracking error can be recovered through convolution of the measured L1-L2 phase difference with the inverse loop transfer function. In this way, corrected L2 observations can be obtained, which are, ideally, free of systematic errors. In practice, the application of the concept suffers from various limitations. First of all, 10 Hz carrier phase observations would be required for the rigorous inversion of a discrete loop with 100 ms update rate rather than the 1 Hz measurements made available in the Swarm mission. Secondly, the finite length of continuous tracking arcs results in edge effects, which limit the quality of the derived phase corrections near the start and end of the track. Despite these limitations, a partial reconstruction of the true L2 phase is possible and helps to improve the overall data quality. Considering periods of high ionospheric activity and narrow loop bandwidth, it is shown that a reduction of carrier phase residuals from 3 mm RMS to 2 mm RMS can be achieved in the reduced dynamic orbit of the Swarm satellites using the proposed correction. Likewise, the artifacts in the gravity field solution around the geomagnetic equator are mitigated, and the quality of the gravity field in the polar regions can be improved. Compared to the down-weighting of erroneous observations used in earlier studies of Swarm gravity field determination, the present correction of looprelated carrier phase errors still shows a larger amplitude of gravity field errors near the equator. However, instead of rejecting or down-weighting observations, the observations can, at least partly, be reconstructed, and a larger number of observations is thus retained. If desired, a combination of a priori loop error corrections with a suitably tuned weighting scheme may be desired for optimum gravity field recovery in future analyses.
Since the systematic errors in the L2 phase observable also affect derived products, for example TEC, a correction should also be considered for other purposes than precise orbit determination. Since the equatorial artifact in GPSonly gravity fields is also known from the GOCE mission, it would be beneficial to know the TL settings for other LEO satellites to investigate if there are also significant artifacts in GPS-only orbits and, if possible, to correct them.