## 2 Methodology and Data
### Weighting of Observation  
Stochastic modeling is presented as VC matrix with variances population only and it depends on the use of the weighting functions. These weighting functions are related to the elevation angle.  
Most of the existing software and studies use the trigonometric functions for weighting observations [58]. In this analysis, the seven weighting functions presented in Table 2 were chosen. First, the model where all observations have the same precision, without a function (hereinafter referred to as Function_1) was chosen. For all other cases, the precision of observations depends on the elevation angle. Function_2 was taken from AC MIT (GAMIT software) [59] and Function_3 from AC National Geodetic Survey (NGS) [58]. Function_4 is adopted from work [60]. Another way to weight observations is to use the threshold value of elevation angle, for example 30\({}^{\circ}\), and to use different functions below and above this angle (Function_5 and Function_6), as was shown in works [11; 32]. The final two functions are basic functions: the simple sine function (Function_7) [61] and square sine function (Function_8) [62]. All these weighting functions were implemented in the PPPH software used.  
PPPH employs the adaptive robust EKF method, which introduces an equivalent weight matrix to compensate for the effect of outliers in observations, and also an adaptive factor to balance the contributions of measurement and estimated parameters [55]. In the [PERSON] filtering method, the functional model and relative dynamic model include the observation noise and process noise with a priori VC matrices, and both are considered to be Gaussian normal distribution (for optimal solutions). In the case of GNSS observation, such an assumption may be incorrect [63; 64]. In our work, the a priori observation VC matrix is a diagonal matrix with variances for code and phase observations for GPS and Galileo observations. The initial uncertainty parameters were determined using the least squares estimation. In this method, the unknown parameters are estimated in the first epoch by the least squares adjustment, and then the estimated parameters and their statistical properties (which may be treated as estimators with Gaussian normal distribution) are utilized as the initial parameters for the next epoch.  
For the calculations, both the pseudorange and carrier-phase observations with set precision of 0.2 m and 0.002 m for raw measurements [18; 65] were used. Because the conventional PPP model uses an ionosphere-free linear combination, it causes a precision of 0.6 m and 0.006 m for the pseudorange and carrier-phase observations of this combination, respectively. Such precision was adopted for the G, E and GE solutions. The Galileo system is still in the phase of building, therefore in the case of the GE1 solution, the precision of 0.3 m and 0.003 m was set for the Galileo code and carrier-phase raw observations, respectively. For this solution, the precision of code and carrier-phase of Galileo observations for ionosphere-free linear combinations are 0.9 m and 0.009 m, respectively. This allowed for an analysis of how the GPS+Galileo solution would be affected by adopting a worse accuracy for Galileo observations.  
The weighting functions used assume the precision of observations is the biggest at the zenith direction and decreases with lower elevation angles. For Function_1, Function_5, Function_6, Function_7 and Function_8, the precision of observation was set as mentioned above. Additionally, in Function_5, the precision of observations increases up to an elevation angle of 30\({}^{\circ}\), and above this elevation angle, precision is the same for all observations. For Function_4, the value was adopted from [35].