## 2 Multi-GNSS ZTD Estimation in Real-Time
For the RT multi-GNSS PPP processing, precise satellite orbits and clocks have to be firstly determined using the observation data from a global GNSS ground tracking network. Similar to the procedure of the IGS ultra-rapid orbits, the RT orbits are predicted (here 6-h prediction) based on the orbits determined in a batch-processing mode by using an orbit integrator. Afterwards, RT satellite clocks are estimated with fixed satellite orbits and station coordinates and updated epoch by epoch due to their short-term fluctuations [19]. The zero mean conditions over the ISB/IFB parameters are also introduced for each system (i.e., BDS and Galileo) and for each GLONASS frequency in the multi-GNSS orbit and clock determination [20]. In the combined data processing, suitable weighting of phase and pseudorange observables among systems should be considered carefully based on their qualities. An elevation-dependent weighting strategy is applied for the combined data processing of four-system observations. The stochastic model of UD pseudorange and carrier phase observations can be described as Equation (3):  
\[ \begin{split}& Cov(i,j)=\left\{\begin{array}{cc}\sigma_{s}^{2}&(i=j) \\ 0&(i\
eq j)\end{array}\right.\\ &\sigma_{s}^{2}=\tau_{s}^{2}\big{(}a_{G}^{2}+b_{G}^{2}\cos^{2}E \big{)},\quad\tau_{G}:\tau_{R}:\tau_{E}:\tau_{C}=1:1.5:2:2\end{split}  \tag{3}\]  
where \(\sigma_{s}\) is the standard deviation of raw measurements (unit: m); \(E\) is the satellite elevation angle (unit: rad); \(a\) and \(b\) are empirical constants. For GPS, \(a\) and \(b\) are generally set to be 0.003 mm for carrier phase and 0.3 m for code observations. For GLONASS, Galileo and BDS, the coefficients are increased by a factor of 1.5, 2.0 and 2.0, respectively, due to the less accurate satellite orbit and clock products compared to GPS. In this study, the empirical weighting factors for different satellite systems are derived from numerous multi-GNSS PPP data processing by using variance component estimation method [21]. Firstly, multi-GNSS PPP is performed by using variance component estimation method, and then the empirical weighting factors for different satellite systems are obtained. Afterwards, the derived weight factors (approximately 1:1.5:2:2) are used in the subsequent real-time PPP processing. The tropospheric slant total delay \(T_{r}^{s}\) consists of the hydrostatic and the wet components, both of which can be expressed by their individual zenith delay and mapping functions,  
\[T_{r}^{s}=M_{h_{r}^{s}}\cdot Z_{h_{r}}+M_{w_{r}^{s}}\cdot Z_{w_{r}} \tag{4}\]  
where \(M_{h_{r}^{s}}\) and \(M_{w_{r}^{s}}\) are the hydrostatic and wet coefficients of the global mapping function (GMF) [22]. The zenith hydrostatic delay (ZHD) \(Z_{h_{r}}\) can be computed by using the [PERSON] model [23] and meteorological data, while the zenith wet delays (ZWD) \(Z_{w_{r}}\), which commonly cannot be corrected well, have to be estimated as unknown parameters. In general, there is no need to estimate tropospheric gradient parameters for the RT ZTD estimation [24]. Additionally, to evaluate the impacts of different treatments of station positions, that is, fixed or estimated as kinematic, on the ZTD estimates, we calculate and make a comparison between the ZTDs derived from the two strategies. A sequential least square filter is employed to estimate unknown parameters in RT processing. All the observations from different GNSS (four systems) are processed together in the common estimator to perform a rigorous multi-GNSS analysis with considering the inter-system and inter-frequency biases [8]. In RT multi-GNSS PPP based atmospheric parameter retrieving, the estimated parameters vector \(\mathbf{X}\) can be expressed as,  
\[\mathbf{X}=(Z_{w_{r}}t_{r}d_{rE}d_{rC}d_{rR_{k}}\mathbf{I}_{r,1}^{s}\overline{ \mathbf{N}}_{r}^{s})^{T} \tag{5}\]  
\[\overline{\mathbf{N}}_{r}^{s}=\mathbf{N}_{r}^{s}+\mathbf{b}_{r}+\mathbf{b}^{s} \tag{6}\]where the receiver clock bias \(t_{r}\) is estimated epoch-wise as white noise.