On the determination of the atmospheric outer scale length of turbulence using GPS phase difference observations: the Seewinkel network

Microwave electromagnetic signals from the Global Navigation Satellite System (GNSS) are affected by their travel through the atmosphere: the troposphere, a non-dispersive medium, has an especial impact on the measurements. The long-term variations of the tropospheric refractive index delay the signals, whereas its random variations correlate with the phase measurements. The correlation structure of residuals from GNSS relative position estimation provides a unique opportunity to study specific properties of the turbulent atmosphere. Prior to such a study, the residuals have to be filtered from unwanted additional effects, such as multipath. In this contribution, we propose to investigate the property of the atmospheric noise by using a new methodology combining the empirical mode decomposition with the Hilbert–Huang transform. The chirurgical “designalling of the noise” aims to filter both the white noise and low-frequency noise to extract only the noise coming from tropospheric turbulence. Further analysis of the power spectrum of phase difference can be performed, including the study of the cut-off frequencies and the two slopes of the power spectrum of phase differences. The obtained values can be compared with theoretical expectations. In this contribution, we use Global Positioning System (GPS) phase observations from the Seewinkel network, specially designed to study the impact of atmospheric turbulence on GPS phase observations. We show that (i) a two-slope power spectrum can be found in the residuals and (ii) that the outer scale length can be taken to a constant value, close to the physically expected one and in relation with the size of the eddies at tropospheric height.


Introduction
The Global Navigation Satellite System (GNSS) refers to constellations of satellites that transmit microwave signals in the L-band. The four GNSS are the European Galileo, the US Global Positioning System (GPS), the Chinese Beidou and the Russian GLONASS system to which can be added local systems such as the QZSS in Japan or NavIc in India (Hoffmann-Wellenhof et al. 2008). The signals transmitted are affected by the Earth's atmosphere; Particularly the signals are delayed and refracted in the troposphere-a non-dispersive medium for microwaves that extends from the ground to approximately 10 km (Stull 1994). Contrary to the ionosphere, its impact on a GNSS signal cannot be mitigated or eliminated by linear combinations (Leick et al. 2015;Teunissen and Montenbruck 2017). The delay in the zenith direction called the Zenith Total Delay has two major components: the Zenith Hydrostatic Delay and the Zenith Wet Delay (ZWD). Whereas the hydrostatic delay can be accurately modelled (see, e.g., Saastamoinen 1973), the ZWD is a major error source in GNSS positioning, particularly for the Up component (Beutler et al. 1987;Dodson et al. 1996). This wet delay is often estimated with the least squares (LS) method as an additional parameter together with the position (Klos et al. 2018), and is an important data source in meteorology (Bevis et al. 1994;Douša and Vaclavovic 2014). It can be used to determine the distribution of atmospheric water vapour for shortterm weather and rainstorm prediction (Champollion et al. 2004;Choy et al. 2013;Li and Deng 2013;Zhao et al. 2018) and for urban applications, such as flooding events (Champollion et al. 2009;Guerova et al. 2016). Long periodic variations of tropospheric wet delays are caused mainly by the steady-state component of the refractive index field.
On the other hand, stochastic variations of the refractive index are induced by tropospheric turbulence that occurs along the signal's path from the transmitter to the receiver (Naudet 1996;Schön and Brunner 2008;Vennebusch and Schön 2011). An accurate knowledge of this second-order effect due to the variations of the refractivity index is of main importance for the higher temporal resolution of water vapour variations. These variations affect microwave phase and phase difference measurements by correlating them (see, e.g. , Hinder 1972;Stotskii and Stotskaya 1992;Buscher et al. 1995 for interferometric observations; Schön and Brunner 2008;Kermarrec and Schön 2014 for GPS;Halsig et al. 2016 for very long baseline interferometry). The observations are no longer stochastically independent: one value recorded at time t will have some similarity with a value taken at t + τ , τ being a time increment and at a different location separated by a spatial distance dr . With a wavelength of approximately 20 cm, GNSS phase measurements are affected by atmospheric turbulence (Tatarskii 1971;Wheelon 2001). Regarding physical considerations, the anisotropic turbulence of the free atmosphere from 1000 m above the Earth's surface is expected to correlate with the phase observations, contrary to the isotropic and well-mixed turbulence of the boundary layer (Myrup 1969;Wheelon 2001).
Residuals from GNSS relative positioning with phase or phase difference observations, e.g., by differencing measurements between receivers to a common satellite, contain precious information about the turbulent atmosphere to assess the second-order effect of the refractive index: a deep analysis of the residuals allows an estimation of some model parameters, such as power law exponents or cut-off frequencies, to be made. The power spectrum density (psd) of the phase difference observations can be physically modelled from theoretical investigations by combining concepts from turbulence theory with electromagnetic wave propagation in a random medium (Ishimaru 2017;Wheelon 2001 and the references therein). Thus, we expect to find a superposition of a fractional Gaussian noise (fGn) with a power law of − 2/3 and a Matérn process with a smoothness of − 4/3 in observation residuals.
Prior to such an estimation from GNSS observations, the geometric part of the observations has to be separated from the stochastic part by adequate processing. A relative positioning scenario, based on so-called double differencing (DD), is often chosen to extract as many systematic effects as possible (Schön and Brunner 2008). Residuals from other positioning scenario, such as Precise Point Positioning or single differencing, need a connection to an ultra-stable external oscillator in order to be interpretable (Schön et al. 2016). Unfortunately, multipath (due to multiple reflections of the signals, see Leick et al. 2015) is responsible for repetitive low-frequency patterns, which bias the power law and cut-off estimation of the tropospheric noise. A mitigation or elimination of such additional effects has to be performed in a preprocessing step to avoid this unwanted effect. We will use the empirical mode decomposition (EMD) in this contribution (Hirrle 2017); as opposed to wavelet decomposition (see, e.g., Satirapod and Rizos 2005), the EMD breaks down the residuals into various components empirically. It allows the fGn (Flandrin et al. 2004) to be extracted with a high reliability because it is independent of predefined basis functions. Thereafter, the correlation parameter of the remaining noise can be estimated and compared with theoretical expectations (see Wheelon 2001, chapter 6).
We will use observations from a network specially designed for estimating turbulent parameters from GPS phase difference residuals to validate the expected theoretical noise model and validate our methodology (Schön and Brunner 2008). The cut-off between the slopes can be accurately estimated using the Whittle maximum likelihood (Sykulski et al. 2019); this allows the computation of the outer scale length of turbulence, which we will interpret based on the separation distance between GPS rays. A further application of this result is the proposal of an adequate stochastic model for GNSS observations to be integrated in the LS adjustment to avoid a problematic overestimation of the precision of GNSS positioning results. We note that the methodology developed in this contribution can be used for other kind of geodetic time series processed using double differenced GNSS phase observations, such as ZWD (Klos et al. 2018).
The remainder of the paper is structured as follows: in a first part, we will briefly present the concept of DD positioning, the second part is devoted to the description of the expected noise, and a third part presents the results

Concept of relative positioning
GNSS satellites transmit microwave signals in the L-band for positioning and timing. The carrier phase signal between a satellite k and a user A can be modelled by where ρ k A denotes the Euclidian distance between the unknown user position and the computed satellite position at the time of signal transmission, c δt A − δt k the synchronization error between the satellite and receiver clock, and T k A , I k A the tropospheric delay and ionospheric signal advance, respectively. Further corrections must be applied for, e.g., antenna phase centre corrections and phase wind up. The wavelength of the signal is denoted by and the ambiguity is N k A , the multipath effect is described by M k A and the observation noise of typical 1-2 mm by ε k A . Interested readers should refer to the basics of GNSS positioning in the textbooks (Hofmann-Wellenhof et al. 2008;Leick et al. 2015;Teunissen and Montenbruck 2017).
Equation (1) underlines the duality between the signal and the nuisance parameter: if the position is known, GNSS signals can be used to study the other terms. Finally, if all terms are corrected or eliminated by combinations of observations, the noise can be studied either with pre-fitted residuals (sometimes called observed minus computed values) or with post-fit residuals after LS adjustment. We will use observed minus computed values in this contribution in order to avoid a potential coupling of the observation residuals due to the adjustment.
Double differences between four carrier phase observations k A , k B , j A , j B are formed to reduce even eliminate error contributions, depending on their similarity-the DD considering two stations A and B and two satellites j and k reads: The elevation of the two satellites at the two stations as well as the atmospheric conditions in the free atmosphere will be nearly identical for medium baseline lengths of up to 10 km (see Fig. 1). Thus, atmospheric delays will be largely reduced, as well as the impact of phase wind und and satellite antenna offsets. Receiver clock errors, satellite clock errors, and hardware delays are eliminated. We mention further that the major part of the first-order ionospheric effect on GPS observation can be removed using the so-called ionosphere-free linear combination of L1 and L2 phase observations (see, e.g., Odijk 2002).
f 2 1 − f 2 2 , where Li , i = 1, 2 indicates measurements in the frequency f 1 = 1575.42 MHz f 2 = 1227.60 MHz. This additional frequency will be used complementary to validate the results.
Finally, the DD time series still contains: • A multipath arising from multiple reflections of the signal on the ground. A far-or near-field multipath will cause mostly low components that can be identified in the residuals (Leick et al. 2015). • Flicker noise (FN) or white noise (WN) comes from electronic devices, cables, receivers or emitters themselves, and • the atmospheric noise (AN) in which we are interested.
In this contribution, we propose to perform a filtering of these additional effects with the EMD and the Page 4 of 16 Kermarrec and Schön Earth, Planets and Space (2020) 72:184 corresponding intrinsic mode functions (IMF). We will follow results from Flandrin et al. (2004) on the fGn, which will be used to model-in its damped versionthe atmospheric or tropospheric noise that affects GNSS phase differences.

Fractional Brownian motion and fractional Gaussian noise
Fractional Brownian motion (fBm) is a self-similar Gaussian process with stationary increments and is a mathematical generalization of the classical random walk and Brownian motion. The term "fractional" is related to the fractional integration and differentiation (Lilly et al. 2017). The fBm has a long-range, strong spatial and temporal correlation. Its increment is called an fGn and was introduced by Mandelbrot and Van Ness (1968) to face the lack of derivative of an fBm. The fGn is a generalization of ordinary WN: its statistical properties are determined by its second-order structure. More precisely, the power law decay of the fGn spectrum is controlled by a single scalar parameter, called the Hurst exponent H , with 0 < H < 1 . For the special case H = 1/2 , the fGn reduces to a WN. For 1/2 < H < 1 , the long-range correlations are said to be persistent, which is the case of interest in this contribution. The power-law exponent β of the spectrum of the process W f ∞1 f β , with f being the frequency, is related to H by H = (β + 1) 2 for an fGn.

Empirical mode decomposition Principle
Huang et al. (1998) developed the EMD to decompose a given time series into a set of components of different time scales called the IMF. Combined with the Huang-Hilbert Transform (hht), this decomposition provides a time-frequency-energy representation of a signal, which is intuitively based on the wave form of the data. It has been used successfully in various domains, from seismic detection to the correction of orbital drift artefacts in satellite data stream, ocean waves or signal analysis, to cite but a few (see Huang and Shen 2005 and the references within). A major field of applications of EMD is to denoise a signal, and various strategies have been proposed to reach that goal. We cite exemplarily regarding geodesy the extraction of WN from GPS co-ordinates time series (Montillet et al. 2013), thresholding-based methods to extract seismic waves (Wang et al. 2012) and mitigation of multipath effects (Wu 2009;Hirrle 2017). The EMD is an iterative method of low computational complexity, in which spurious harmonics are less likely to appear: no linear property is imposed on a possibly nonlinear system. Its empirical approach is independent of mother wavelets-as in wavelet analysis-or of some trigonometric functions-as in Fourier transform. The EMD when used on non-stationary signals, allows a direct extraction of the local energy associated with the IMF. Its main drawback is linked to the non-uniqueness of the decomposition; the latter is related to a particular interpolation scheme for envelope extraction, the stopping criterion (usually 1e −6 ) used in the sifting process and the number of IMF chosen.
The EMD algorithm is summarized in Flandrin et al. (2004) as the following loop: • identification of the extrema of the signal. • interpolation between minima (maxima) to get an envelope of the signal. • computation of the average of the envelope.
• extraction of the detail as the difference between the envelope and the original signal. • iteration on the residuals.
An IMF should also have the following properties: (i) zero mean, (ii) all the maxima and all the minima of one IMF are correspondingly positive and negative, and (iii) narrowband. Each IMF has fewer extrema than all the lower order IMFs, leading to decreasing oscillations as the IMF order increases.
The original temporal signal S(t) is finally a nonweighted combination of all the IMF obtained.
where d i (t) is the ith IMF and m K (t) the residuals after K decomposition.
A variance σ 2 IMF,i , corresponding to its mean energy, is associated to each IMF. Under "energy", we understand the empirical variance estimate of a given IMF.

Intrinsic mode functions of a fractional Gaussian noise
The IMF of an fGn and their corresponding energies are studied in detail in Flandrin et al. (2004) and Wu and Huang (2004) for the specific WN case ( H = 1/2 ). The authors estimated the psd of each mode using dedicated simulations and showed that the EMD is acting as a dyadic filter bank of constant-Q band pass filters on the fGn.
Apart from the first noise-only IMF, which carries most of the WN component and, thus, the main energy, the variance of the other IMF d i (t), i = 2 . . . K decreases linearly in a semi-log diagram regarding the IMF index.
Decomposition of the fGn: slope of the IMF variance To illustrate this property, we simulated three fGn, each of length 60,000, corresponding to a Hurst parameter of 0.5 (WN), 1 (FN, red line) and 0.83 (atmospheric fGn, AN, see "The atmospheric noise from phase difference" sec- Page 5 of 16 Kermarrec and Schön Earth, Planets and Space (2020) 72:184 tion or Wheelon 2001). Please note that we will use these three abbreviations for the names of the different noises hereafter and for the sake of readability. We used the Matlab function ffgn to generate the time series (Stoev 2020). The latter are depicted in Fig. 2 (left), together with the IMF log 2 variance ( Fig. 2 right) estimated with the Matlab function hht and emd using seven decomposition steps. The theoretical expectation that the IMF variance should decrease exponentially as a function of the IMF index is confirmed. The corresponding negative slope κ H , starting from the second IMF, gives an estimation Ĥ of the Hurst parameter: In the following, we will use the following notation: IMF i is the ith IMF, i the index of the ith IMF, i.e. i = [1, . . . , K ] , so that the log 2 variance of the ith IMF can be simply expressed as: C being a constant. Figure 2 (right) highlights further that: • The slope of a noise combination consisting of AN and WN (magenta line) is between the one of a single WN (yellow line) and AN (red line) from IMF 2 . • The slope for an FN is nearly 0 from IMF 2 and is − 1 for a WN. Thus, when the total noise is a combination of power law noises, the slope of the log 2 vari- ance gives an indication of the dominant source, which may be FN, WN or AN. • An AN "diluted" in a combination of WN and FN of equal variances (green line) leads to a change in the slope regarding a simple AN (red line) for a finite sample size. This change will depend on the respective variances of the different noises. The slope of the green line in our example is affected by the WN component up to IMF 5 . • The log 2 variance of IMF 1 always behaves differently than the one of higher index and should not be integrated in the slope estimation. This is coherent with the empirical results of Flandrin et al. (2004Flandrin et al. ( , 2005.

Proposal to extract atmospheric noise First step: elimination of multipath/flicker noise signal
Any discrepancies from the linear relationship between the mean energy of the IMF and their index can be used to denoise a signal (see, e.g., or Wu and Huang 2004;Kopsinis and McLauglin 2009), as highlighted in Fig. 2 (right). In this contribution, however, we wish to analyse the noise and not the signal. We will use the terminology "designalling the noise" intentionally where the signal comes here from a multipath signature, additional noise, such as FN, or noise from functional misspecifications. Such a model was chosen to compute the position. The first step of the procedure to extract the AN consists of estimating the slope κ H from the mean energies of IMF i , i ∈ [2, 3, 4] , which is performed by regression. Based on the linear slope parameter κ H and C obtained (see Eq. (5)), we compute the standard error for Page 6 of 16 Kermarrec and Schön Earth, Planets and Space (2020) 72:184 prediction , for which κ H ± � corresponds to a roughly 68% prediction interval. The strategy to filter the signal from the noise by excluding some IMF was firstly develop for WN. In this contribution, we extend it to an fGn, following the proposed methodology in Wu and Huang (2004) for WN. Consequently: • We iteratively control whether log 2 σ 2 IMF,i = log 2 σ 2 IMF,i ± 3�, i = 5 . . . K to decide if the ith IMF should be further considered or not. If not, the corresponding IMFs are no longer dominated by an fGn but rather by multipath, functional misspecifications or eventually noise with a higher power law than the AN (see Fig. 2, green line), i.e. what we called the "signal". • The filtered S filt (t) is finally computed as: As suggested in Fig. 2 (right), the slope estimated may be affected when the noise is an unknown combination of different power law noises, as exemplarily, AN + WN + FN. Following results from device design and analysis (Bahreyni 2008), we consider that the FN is to be found at low frequencies, i.e. in the higher IMF, and, thus, eliminated by the previous filtering. Such results are coherent with the previous comment made in Fig. 2. It remains, thus, to filter the WN from the remaining fGn S filt (t) in order to extract the AN only, which we perform using the EMD.

White noise filtering
From the last step, the high-pass filtering led to . We assume that this remaining noise can be decomposed into two components: is the WN and S AN (t) the AN. Unfortunately, the WN component affects the power spectral density (psd) of S filt (t) strongly at high frequencies by driving its slope toward 0: an interpretation of S filt (t) as an AN is made meaningless without a prior whitening. In a first step, we start the filtering using the empirical maximization algorithm developed in Kargoll et al. (2018) in order to determine the ratio R WN of the standard deviation of the WN to the one of S filt (t) . We further assume that the main part of the WN is present in the first IMF of the EMD decomposition, which contains most of the variance of S filt (t) (Flandrin et al. 2004;Wu and Huang 2004;Taebi 2017). Consequently, . The AN part is, thus, finally written as: This component will be further analysed based on the physical consideration developed in "The atmospheric noise from phase difference" section.

Summary of the methodology to extract the atmospheric fractional Gaussian noise
The final procedure to extract the AN is summarized in a flowchart form in Fig. 3. It involves three steps: 1. Computation and preprocessing of DD residuals (mean and trend extraction). 2. Elimination of the non-noise component (mainly multipath and low frequency FN).
3. Filtering of the WN component by downweighting the first IMF.

Further comments
Following Flandrin et al. (2004Flandrin et al. ( , 2005 and without loss of generality, we will fix K = 7 . We make use of the Matlab built-in function emd. The interpolation method chosen is pchip, which is dedicated to quick varying signals.

The atmospheric noise from phase difference
Once the AN has been extracted from the DD residuals, its correlation structure and phase spectrum can be analysed based on a physically plausible consideration. Consequently, we follow Wheelon (2001, chapter 6, p. 304-305). In a first subsection, we will briefly present the assumptions made, followed by a description of the two slopes noise model and its dependencies.
Page 7 of 16 Kermarrec and Schön Earth, Planets and Space (2020) 72:184 The Matérn process or a damped fractional Brownian motion to model atmospheric noise

Atmospheric temporal correlations
Since turbulence theory is usually concerned with the statistical moments of the velocity or its fluctuations. Turbulent variations of the refractive index in the atmosphere can be modelled using the concept of eddies or swirl of motions (Stull 1994). These energy structures develop with wind shear, which is most often caused by microbursts from thunderstorms, temperature inversions, and surface obstructions. Eddies may have horizontal or vertical circulation. In the Earth's atmosphere, small eddies occur in the daytime when surface air is warmed by contact with the ground and then rises. On a larger scale, the high-and low-pressure centres that move over the Earth may be considered eddies, i.e. they are anomalous to the general circulation. Consequently, the turbulent process is described using an energy cascade in which eddies subdivide into progressively smaller eddies until they are dissipated by viscosity or heat. The cascade starts at the outer scale length wavenumber and ends when eddies have the size of the inner scale length. Phase measurements are mainly influenced by large eddies, so that the main region of interest in the cascade is the energyinput region where the atmospheric wind field becomes unstable and creates large, two-dimensional energetic structures (Wheelon 2001). Figure 4 proposes a schematic representation of the troposphere where eddies of different size and energy coexist. Eddies in the boundary layer are more likely to be small and isotropic with a high reorganization rate: the so-called three-dimensional turbulence. On the contrary, anisotropic eddies, elongated mainly in the horizontal direction, will be present in the loosely called free troposphere, from the tropospheric height at approximately 1000 m altitude (Stull 1994).
Using a geometric optics models and having in mind a positioning with DD processing, a ray coming from a GNSS satellite is expected to encounter the same elongated eddies as another ray coming from the reference satellite. This assumption holds true as long as the distance between the rays do not exceed the size of the largest eddies that affect the phase measurements. Figure 5 explains this concept schematically, highlighting that the rays that are spatially close (bold dotted and dotted black for the non-reference "green" satellite) will be temporally correlated. We call that effect "temporal correlations with spatial components".

Modelization of the power spectrum for phase differences
Many natural noises have approximate isotropic power spectra of the power-law form. From dimensional analysis, Kolmogorov (1941) has shown the power spectrum of passive scalars, such as temperature or refractive index can be written as: Fig. 3 The methodology to extract an AN in three steps Page 8 of 16 Kermarrec and Schön Earth, Planets and Space (2020) 72:184 where C 2 n is the structure constant of the atmosphere and κ = [ κ x κ y κ z ] T , the vector of three-dimensional wavenumbers. This model is valid in the inertial range for homogeneous and isotropic turbulence, where 2π L 0 ≤ κ ≤ 2π l 0 with L 0 , l 0 being the outer and inner scale length of turbulence that bounds the inertial range, respectively.
In order to preserve the behaviour of the spectrum in the inertial range, the empirical Von Karman model is often preferred to the Kolmogorov description. It reads: , with κ 0 being the wavenumber corresponding to a cutoff frequency that allows finite expressions for the first two moments of the spectrum. As mentioned previously, anisotropic eddies are more likely to correlate with the GNSS signals. Stretched rectangular co-ordinates can be used to account for the anisotropy and the spectrum is expressed as: where a, b and c are stretching factors describing the elongation of the eddies in the three dimensions. Usually, a = b describes the horizontal elongation of the eddies, whereas c is related to the vertical stretching. The residuals of the DD processing of GNSS phase observations are time-dependent. Consequently, only the temporal correlations can be assessed: the spatial phase covariance between two GNSS rays from two satellites (10) � n (κ) = abc� n a 2 κ 2 x + b 2 κ 2 y + c 2 κ 2 z , Boundary layer: Height < 1000 m Free atmosphere: elongated eddies, Height > 1000 m can be transformed into its temporal counterpart by means of the co-ordinate transformation r = r ′ − v eddies τ using Taylor's frozen hypothesis (Taylor 1938), r ′ and r being the position and the transformed position, respectively. With this approximation, eddies are considered to be frozen during the measurement, only moved by the mean wind v eddies , which is further assumed to blow parallel to the baseline. It can be replaced by a mean velocity of v eddies = 8 ms −1 , a good approximation at the top of the boundary layer (Stull 1994). This value corresponds to the geostrophic wind, which is a theoretical wind balanced by the Coriolis and pressure gradient forces. Under the approximation mentioned previously, we call ρ ray , i.e. the Euclidian distance between 2 GNSS rays from the same satellite, computed at a given height. Following Wheelon (2001, chapter 6) or Tatarskii (1971, chapter 4) and using the previous assumptions, the power spectrum of phase difference = k − j between the phases of two satellites k and j can be derived from the single-path power spectra W � (ω) as where ω is the angular frequency. The factor sin 2 ωρ ray 2v eddies corresponds to the interference pattern of the receiver (Ishimaru 2017). By integrating along the lines-of-sight, W � (ω) is expressed as the product of a factor that accounts for the elevation of the satellite and a damped power law spectrum, i.e.
where α = κ 0 v eddies a = 2π v eddies aL 0 , El is the elevation of the satellite, H trop the tropospheric deep or height taken to 1000 m (Stull 2009), a = b and c the horizontal and vertical stretched parameters for the outer scale length, respectively, and C 2 n is the structure constant, which varies with height but is here considered as a constant. k is the electromagnetic wavenumber. W � (ω) corresponds to the power spectrum of a Matérn process, i.e. a damped fBm (Lilly et al. 2017). The cut-off frequency α = 2π v eddies aL 0 is related to the wind velocity in the free atmosphere and the horizontal outer scale length of turbulence aL 0 , which is known to be several kilometres along most of the path travelled by the GPS signals from the satellites to the receiver. We may conjecture from Fig. 3 that the cut-off frequency could correspond to the maximum size of the eddies that will correlate with the GPS signals between two epochs. The threshold angular frequency α = 2π v eddies aL 0 is less than 0.01 Hz for (11) W �� (ω) = 4 sin 2 ωρ ray 2v eddies W � (ω), W � (ω) = 2.192H trop k 2 C 2 n ca −5 / 3 v 5 / 3 eddies sin 2 (El) 1 ω 2 + α 2 5 / 3 , v eddies = 8 ms −1 . When the time series of residuals are long enough to allow for a trustworthy estimation can only be estimated, i.e. at least 5000 epochs a 1 s. Following Wheelon (2001) and from Eq. (12), the spectrum of GNSS phase difference should behave as ω − 2 / 3 below a corner frequency of ω c = 2v eddies ρ ray , i.e. for sin 2 ωρ ray 2v eddies ∼ ω ω c 2 . Consequently, we expect to have two slopes in the psd of GPS phase difference residuals: Using this decomposition in three domains, we consider that the Matérn process corresponding to W will start at a high frequency as an fGn and will switch to a Matérn model at a lower frequency, i.e. a saturation.
Note on double differences we assume that the power spectrum of phase double differences in this contribution is the same as the power spectrum of phase differences. By making this approximation, we consider a "mean" station between A and B and take the difference between the phase from the same satellite k (or j) at two stations, i.e. �(t) = � k B − � k A (t) , as the phase measurement. The value ρ ray is not the distance between the two stations. We consider that the latter is not playing a role in the correlation structure, which we justify by the relatively long baseline length chosen in our case study regarding the value of the outer scale length of turbulence expected.
Note on ionospheric correlations ionospheric correlations due to plasma irregularities can be modelled with a similar power law spectrum as tropospheric turbulence (Rino 1979). Using GPS observations from Swarm satellites, Kermarrec et al. (2018) demonstrates that a similar model to that for tropospheric correlations can be used for ionospheric irregularities. However, the outer scale length is expected to be up to five times larger than in the free atmosphere and cannot be detected with the relatively short residuals samples of DD adjustment. That is the reason we allow ourselves to loosely use the terminology "atmospheric correlations" to describe correlations due to the turbulent variations of refractive index in the troposphere only. For the sake of completeness and to validate this approximation, we will use observations from the ionospheric-free linear combination in our real data analysis (Leick et al. 2015).
Page 10 of 16 Kermarrec and Schön Earth, Planets and Space (2020) 72:184 How to estimate the noise parameters

The debiased Whittle maximum likelihood
Dealing with real observations, the noise parameters of the damped fBm will have to be estimated from the preprocessed DD residuals. Three parameters of the correlation function have to be estimated from Eq. (12): the variance, the range and the smoothness. This estimation can be done with the maximum likelihood estimator: (i) as in the so-called Hector software (Bos et al. 2008) with an adapted version of the Matérn model-the Gauss Markov model-or (ii) with the powerful Whittle maximum likelihood estimator. This method leads to a fast and nearly unbiased estimation of the noise parameter for large samples (Sykulski et al. 2019). It further produces estimates that are asymptotically efficient and avoids biases coming from visual slope determination. Usual ML strategies involve matrix inversion: they are avoided when the autocovariance matrix is approximated in the frequency domain. Unfortunately, the so-called "Whittle likelihood" uses the periodogram of the observations, which may be a biased measure for finite samples due to blurring and aliasing effects (Percival and Walden 1993). Sykulski et al. (2019) proposed an alternative based on quasi-likelihood called the "de-biased Whittle likelihood", which is particularly interesting for discretetime processes of finite length, such as the residuals of the LS adjustment used in this contribution. The Whittle likelihood in its discretized form is given by , with the set of discrete Fourier frequencies, f (ω, H ) the continuous-time process spectral density and I(ω) the periodogram I(ω)∞ N j=1 z j e −ijω 2 . The Whittle estimator assumes a priori that the power spectrum of the underlying process of the dataset is known and should only be used if a time series is consistent with the process under consideration. In a first approximation, the visual inspection of the psd validates the presence of typical Matérn damping at low frequencies (see Lilly et al. 2017) or not.
In this contribution, we use the Whittle maximum likelihood estimator as implemented in Matlab by Lilly (2020).

Parameters to be estimated
From Eq. (12), we have to force for the estimation of a Matérn noise (see Lilly et al. 2017) corresponding to a smoothness of 4 3 and a cut-off frequency ω c = 2v eddies ρ ray . Furthermore, a second Matérn noise with a smoothness close to 1 3 should be estimated from the residuals with a saturation at α = 2π v eddies aL 0 . This value is expected to be close to a constant, as it depends mostly on atmospheric parameters, or at least on the maximum size of the eddies that will correlate with the GNSS measurements. Unfortunately, the Matérn process is only defined for a smoothness higher than 1 2 , so that the estimator is constrained for a damped FN estimation.

Real data analysis
In this section, we validate the methodology using GPS observations from the Seewinkel network (Schön and Brunner 2008;Kermarrec and Schön 2014). We note that our methodology is not dependent on the GNSS system and can be extended for Galileo, Beidou, GLONASS or further satellites systems, provided that phase observations are used and processed using a relative positioning.
In this case study, we will show that an estimation of the cut-off frequency from the DD residuals is (i) feasible and (ii) provides physically plausible estimation. We further aim to analyse its dependencies.
The network will be presented briefly in a first subsection. A second part explains didactically the methodology and the computation of the outer scale length and cut-off frequency for one time series of residuals, which is extended to all visible satellites in a third subsection.

Observations
We use GPS phase observations from the Seewinkel network, which was specially designed to study the impact of temporal correlations due to turbulence (Schön and Brunner 2008). This network is located in a flat terrain in eastern Austria (Burgenland,47.7 N,17.0 E, 160 m ellipsoidal height) and consists of six aligned stations P0, P1, P2, P4, P8 and P16 with separations of approximately 1, 2, 4, 8 and 16 km, respectively, measured on 15 April 2003 during 8 h (5:45-13:45 GPS time) using identical equipment, a 1 Hz data rate and a cut-off-angle of 3° (see Fig. 6). We use the stations up to P8 so that the assumptions that DD are eliminating ionospheric refraction still holds. In this contribution, we use L1 observations (1575.42 MHz signal) without loss of generality and the ionospheric-free combination obtained from L1 and L2 (1227.60 MHz signal) additionally (Leick et al. 2015).
The GPS processing was done with the Bernese GPS Software Version 5.0 (Dach et al. 2009). In a first step, the coordinates of the stations P0 to P16 as well as relative troposphere parameters were precisely determined as a relative network w.r.t P0 using the ionosphere-free linearcombination L0. The ambiguities were fixed using the Fig. 6 The alignment of the stations of the Seewinkel network Page 11 of 16 Kermarrec and Schön Earth, Planets and Space (2020) 72:184 SIGMA strategy in a widelane-narrow-lane run. The coordinates of P0 were held fixed to a previously determined position w.r.t the GNSS station in Graz. In these computations, the whole available data material of ca. 8 h is used with a data rate of 30 s. Precise final orbit and clocks products from the International GNSS Service (IGS, Johnston et al. 2017) were introduced as well as antenna correction. In a second step, the coordinates of the stations were held fixed and the undifferenced OMC for different frequencies (L1, L2, L0) were computed sequentially for each station with the Bernese software for a 1 s data rate. Here, the 5-s CODE satellite clocks, final orbits and global ionospheric model were used (Dach et al. 2020). The receiver time scale was synchronized to GPS time by code observations and all typical corrections were applied, like antenna phase centre correction for user and satellite antennas from the IGS. Using a specially implemented writing option, the raw carrier phase observations as well as all corrections computed by Bernese were stored for 2 h data segments. This length allows for a reliable estimation of the correlation parameter.
In a last step-in own MATLAB code-the double differences were formed following Eq. (2), and the ambiguities solved by rounding the OMC to the nearest integer.
For the sake of shortness, we only present the results corresponding to P1 and P8, for which the North, East and Up (N, E, U) components are estimated. The skyplot corresponding to the observation time span is presented in Fig. 7. PRN6 was taken as the reference satellite for the DD computation (PRN or pseudorandom noise can be considered as the satellite number). An a priori standard deviation of 1 mm was assumed for the L1 carrier phase measurements.
Although the GPS observations may not correspond to the actual state of the art, the network is unique as the stations are ideally aligned. It allows to study the cut-off frequency predicted from atmospheric turbulence theory, which its designer Prof. F.K. Brunner was searching for. Further contribution will focus on more actual data set using, e.g., single differences.

Didactical example Preprocessing of the DD residuals
We use the DD between PRN5 and PRN6 from baseline P0P2 for our didactical example. This first batch of observations consists of 5050 epochs during which the satellites were visible. The elevation of PRN5 is between 40° and 5°. The corresponding DD residuals are plotted in Fig. 8 (left). From epoch 3000 s, we recognize the typical sinusoidal envelope corresponding to a strong far field multipath effect due to the decreasing elevation of the satellite: it has to be eliminated prior to the AN analysis.

Elimination of signal components
After preprocessing the DD, the second step consists of extracting the noise from the multipath signal. Consequently, we firstly compute the EMD of the preprocessed DD residuals (see Fig. 8, right). The visual inspection confirms that the first four IMFs are more likely to contain the high-frequency noise of the residuals; the latter may be white, coloured or a combination of both. On the contrary, higher modes should correspond to the multipath envelope or unmodelled periodic component.
After having decomposed the noise into IMF, we perform the hht transformation and compute the log 2 of the mean of the energy for each IMF. We estimate a slope of − 0.95 for the line between the log 2 variance between indices 2-4 by regression (Fig. 9, left, red line). The corresponding domains of exclusion corresponding to 3Δ are additionally depicted as dotted blue lines: They allow us to include the IMF 5 in the noise, whereas IMF 6 and 7 have to be eliminated prior to further analysis. The slope is recomputed by adding IMF 5 and reaches − 0.90, highlighting a high contamination with WN. However, the high negative slope does not speak in favour of a high contamination with FN at high frequencies.

Filtering the WN
In a second step, the multipath-free time series is further filtered by means of an adequate recombination of the EMD. A coefficient of R WN = 0.35 was estimated. The time series obtained S AN (t) is plotted in Fig. 9 (right, up) together with its psd estimated using the pwelch function of Matlab (Fig. 9, right, bottom).  Kermarrec and Schön Earth, Planets and Space (2020) 72:184 We identify a psd corresponding to a Matérn process visually, by the typical damping at low frequency. Two slopes can be clearly identified, as well as their corresponding cut-off frequencies.
We use the debiased Whittle maximum likelihood to estimate these parameters. We found an outer scale length L 0 = 3900 m and a separation distance ρ ray = 5.1 m, corresponding to the transition − 2/3 and − 8/3.

Interpretation of the results
We note that the outer scale length of turbulence corresponds to a plausible size of the eddies at the tropospheric height. At that height, we guessed, from physical consideration on electromagnetic wave propagation in random medium, that turbulence affects the GNSS signals. We also estimated the mean Euclidian distance between the GPS rays of PRN 5 and the reference satellite (PRN 6) at that height during the batch of 5050 epochs under consideration. The mean separation distance Page 13 of 16 Kermarrec and Schön Earth, Planets and Space (2020) 72:184 reaches a value of 3310 m, which is coherent with the estimated cut-off length L 0 , i.e. the maximum size of the eddies that affect GPS phase differences. The cut-off separation distance at high frequency is slightly more difficult to interpret but could be linked, following Wheelon (2001 chapter 6), with the Euclidian distance between the rays from PRN 5 and 6 near the ground. This value is approximately ρ eff = 4.9 m at a height of 2.5 m above the ground, a value close to the cutoff distance estimated. Clearly, this is only an approximation as in that case v eddies should not be approximated by the geostrophic wind velocity, leading to an even smaller value of ρ ray .
Note on the estimation error an error of �α = 0.01 for log 10 (α) = −1 translates to |�L 0 | = 2π v eddies α 2 �α ≈ 50 m, which is less than 5% of the nominal value, considering L 0 around 4000 m. We consider this uncertainty as acceptable, regarding the unknown and relatively high value of the outer scale length. An error of �α = 0.01 for log 10 (α) = 0 translates to | L 0 | ≈ 0.5 m, i.e. ρ ray is slightly more difficult to estimate with a high accuracy.

Baselines P0P2 and P0P8 Results
In this section, we propose to extend the analysis of the preprocessed residuals for all visible satellites using observations from the first and second batch of 2 h of observations for two stations P2 and P8. These two stations were chosen to compare the impact of increasing distance on the parameters estimated, i.e. to validate the assumption made in "Modelization of the power spectrum for phase differences" section about the nondependency of the results regarding the baseline length. Other results are not presented due to their similarities and for the sake of briefness, but can be made available on demand.
We follow the two-step methodology detailed in the previous section: (i) extracting the multipath component and (ii) eliminating the WN component. The results are presented in Table 1. We add the results of P0P2 for the ionospheric-free linear combination in Table 2.

Comments
Outer scale length of turbulence L 0 Table 1 highlights that the outer scale length of turbulence L 0 can be estimated between 3500 and 3900 m from the GPS double difference phase power spectrum. This value corresponds approximately to the median value over all satellites and epochs of the separation distance, which reaches 4050 m for the first batch of observations and 4120 m for the second one, computed at the tropospheric height of 1000 m. Thus, a plausible interpretation is that L 0 may correspond to the maximum size of eddies that act as correlating the phase difference. The low-frequency cut-off parameter is dependent neither on the elevation of the satellites nor on the baseline length, which is physically followable.
The ionospheric-free linear combination (Table 2) does not show any discrepancy with the previous results: the estimated mean value of L 0 is close to the  Kermarrec and Schön Earth, Planets and Space (2020) 72:184 value found without eliminating the impact of the ionosphere. This confirms that the outer scale length in the ionosphere is probably much higher than in the troposphere, so that it cannot be detected in the relatively short time series. Moreover, its impact on the phase difference, mainly due to plasma irregularities (Rino 1979), is much smaller in magnitude than the one of the tropospheric turbulence.
The separation distance or cut-off at high frequencies The cut-off ρ ray reaches values between 4 and 8 m for all satellites, and its mean is close to 5 m. It could be related to the separation distance between two rays computed near the ground, at a height of 1-2 m. It has some similarities with values found for optical turbulence (Coulman and Vernin 1991). However, we propose not to over-interpret this value, which is difficult to estimate with a high reliability since it does not correspond to a saturation of the power spectrum: the slope transition is not straightforward to identify. A bound within the expected range of frequencies for the estimation with Whittle maximum likelihood is necessary. Visually, ρ ray can be considered as nearly constant for all baseline lengths and satellites. It is not affected by the elimination of ionospheric effects (Table 2), as the latter should mainly be detected at low frequencies.
Although the values of Table 2 are slightly  higher than for Table 1 with ionospheric impact, they should be weighted by the fact that the ionosphere-free linear combination contains a strong WN component. The latter cannot be fully eliminated and impacts the determination of the cut-off at high frequencies.
The non-dependency with the baseline length is coherent with our assumption that a mean phase can be considered for medium baseline lengths (up to 30 km) to compute the phase spectrum of the double differences. The low dependency with the elevation of the satellite is justified by the small batch lengths, making it difficult if not impossible to calculate an accurate temporal or, at least, batch-wise determination of this parameter. Further analysis and experiments are planned with highquality observations from single differences to validate a potential elevation dependency, as well as dependency with the wind vector.

Conclusion
A knowledge of the stochastic properties of GNSS phase difference observations is mandatory to get a reliable estimation of the position and for further processing implying test statistics. Any misspecification can lead to over-or underestimation of the precision. The study of the noise can be performed by analysing the residuals from DD processing of GNSS phase observations. In addition to WN or/and FN, it is expected that atmospheric noise will correlate with the observations for medium baseline lengths. These typical correlations come from the path of the GNSS signals through the atmosphere. They result from plasma irregularities in the ionosphere and turbulent structures called eddies in the troposphere. The latter source is considered to be stronger for GNSS phase observations and occurs in the loosely called free atmosphere. There, the reorganization of the eddies is less intense than near the ground, remaining, however, linked to strong values of the structure constant and a non-negligible energy. Using concepts from the turbulence theory, the psd of phase differences for microwaves in the turbulent troposphere can be modelled as a Matérn model modulated by a sinus term. Consequently, the spectrum of the phase difference residuals should present two slopes: • − 8/3 at high frequencies followed by • a saturating − 2/3 power law noise. The cut-off frequency of saturation corresponds to that of a damped Matérn process. This value can be related to the outer scale length of turbulence, i.e. the maximum size of the eddies acting as correlating the phase measurements, from physical consideration of the anisotropic turbulence in the free atmosphere.
These parameters can be estimated using the debiased Whittle maximum likelihood, which has the main advantage of being less biased for short samples than the normal maximum likelihood estimator. Estimating the parameters, as, e.g., the outer scale length leads to a more reliable results than finding slope intersections visually.
Prior to the estimation of the noise parameters, the residuals have to be filtered from WN and additional low frequency effects, such as multipath. In this contribution, we developed a strategy to extract the tropospheric noise from phase DD residuals by revisiting the EMD in an alternative way, which could be called "designalling the noise". We successfully applied our methodology to GPS observations from medium baselines from a network called "Seewinkel", specifically designed for the analysis of the impact of turbulence on GPS phase observations from DD phase residuals. We could estimate the threshold frequencies from a saturating two-slopes power spectrum. On the selected day and selected area, there were dynamic conditions on the periphery of upper-air anticyclone to the west. Given those conditions, the obtained results for outer scale turbulence between 3.5 and 4 km are plausible. We found furthermore that this value was related to the separation distance between GNSS rays, computed at the tropospheric height, where the main impact of the turbulence on GNSS observations is expected. Our findings validate that: • the extraction of the tropospheric noise from GPS phase residuals can be performed with a high reliability; the ionospheric-free linear combination confirmed the result, • it is possible to derive physically plausible turbulence parameters-such as the outer scale length of turbulence-that can be linked with meteorological conditions, • the stochastic model for GNSS phase measurements can be modelled, i.e. without empirical iterative computation, thus avoiding the overestimation of the precision.
We further computed the cut-off at high frequencies. The value obtained was found to be nearly constant for all satellites and baseline lengths and close to 5 m, a value that is probably related to the separation distance computed near the ground.
Besides being useful in meteorology, GNSS phase observations can be used to analyse the turbulence in the atmosphere, with the restriction that the parameters may not be as accurate as with dedicated sensors.
Further analysis is planned to confirm our finding by analysing data sets from high accurate single difference observations. We expect to find a dependency of the cut-off frequency on the elevation of the satellites by a batch-wise estimation of the parameters.