## 2 Multi-GNSS ZTD Estimation in Real-Time  
In the RT multi-GNSS PPP processing, the linearized equations of undifferenced (UD) carrier phase and pseudo-range observations can be expressed as follows,  
\[\left\{ \begin{array}{l}P_{r,j}^{\rm G}=-\mathbf{u}_{r}^{\rm G}\cdot x+{t_{r}} \,^{-\rm G}-t^{\rm G}+\kappa_{jG}\cdot I_{r,1}^{\rm G}+T_{r}^{\rm G}+\lambda_{ jG}(N_{r,j}^{\rm G}+b_{rG,j}-b_{j}{}^{\rm G})+\epsilon_{r,j}^{\rm G}\\ {I_{r,j}^{\rm R_{k}}}=-\mathbf{u}_{r}^{\rm R}\cdot x+{t_{r}}\,^{-\rm R}-t^{\rm R }+\kappa_{jR_{k}}\cdot I_{r,1}^{\rm R}+T_{r}^{\rm R}+c\cdot d_{rR_{k}}+e_{r,j}^ {\rm R}\\ {P_{r,j}^{\rm E}}=-\mathbf{u}_{r}^{\rm E}\cdot x+{t_{r}}\,^{-\rm E}-t^{\rm E}+ \kappa_{jE}\cdot I_{r,1}^{\rm E}+T_{r}^{\rm E}+c\cdot d_{rE}+e_{r,j}^{\rm E}\\ {P_{r,j}^{\rm C}}=-\mathbf{u}_{r}^{\rm C}\cdot x+{t_{r}}\,^{-\rm C}-t^{\rm C}+ \kappa_{jG}\cdot I_{r,1}^{\rm C}+T_{r}^{\rm C}+c\cdot d_{rC}+e_{r,j}^{\rm C} \end{array} \right. \tag{1}\]  
\[\left\{ \begin{array}{l}L_{r,j}^{\rm G}=-\mathbf{u}_{r}^{\rm G}\cdot x+{t_{r} }\,^{-\rm G}-t^{\rm G}-\kappa_{jG}\cdot I_{r,1}^{\rm G}+T_{r}^{\rm G}+\lambda_{ jG}(N_{r,j}^{\rm G}+b_{rG,j}-b_{j}{}^{\rm G})+\epsilon_{r,j}^{\rm G}\\ L_{r,j}^{\rm R_{k}}=-\mathbf{u}_{r}^{\rm R}\cdot x+{t_{r}}\,^{-\rm R}-t^{\rm R }-\kappa_{jR_{k}}\cdot I_{r,1}^{\rm R}+T_{r}^{\rm F}+\lambda_{jR_{k}}(N_{r,j}^ {\rm F}+b_{rR,j}-b_{j}{}^{\rm G})+\epsilon_{r,j}^{\rm R}\\ L_{r,j}^{\rm E}=-\mathbf{u}_{r}^{\rm E}\cdot x+{t_{r}}\,^{-\rm E}-t^{\rm E}- \kappa_{jE}\cdot I_{r,1}^{\rm E}+T_{r}^{\rm E}+\lambda_{jE}(N_{r,j}^{\rm E}+b_ {rE,j}-b_{j}{}^{\rm E})+\epsilon_{r,j}^{\rm E}\\ L_{r,j}^{\rm C}=-\mathbf{u}_{r}^{\rm C}\cdot x+{t_{r}}\,^{-\rm C}-t^{\rm C}- \kappa_{jC}\cdot I_{r,1}^{\rm C}+T_{r}^{\rm C}+\lambda_{jC}(N_{r,j}^{\rm C}+b_ {rC,j}-b_{j}{}^{\rm C})+\epsilon_{r,j}^{\rm C}\end{array} \right. \tag{2}\]  
where the indices \(G\), \(R\), \(E\) and \(C\) refer to the GPS, GLONASS, Galileo, and BDS satellites, respectively; \(R_{k}\) expresses the GLONASS satellite with frequency factor \(k\); \(s\), \(r\), and \(j\) represent satellite, receiver, and frequency, respectively; \(P\) and \(L\) denote \"observed minus computed\" pseudorage and phase observables, respectively; \(\mathbf{u}_{r}^{\rm s}\) is the unit vector of the direction from receiver to satellite; \(x\) is the positioning error of station coordinates; \({}^{s}\) and \(t^{\rm s}\) refer to the RT satellite orbit errors and satellite clock errors, respectively; \(t_{r}\) denotes the clock errors of receiver; \(d_{r,s}\) is the code biases; \(N_{r,j}^{s}\) is the integer ambiguity; \(b_{r,j}\) and \(b_{j}{}^{s}\) are the uncalibrated phase delays; \(\lambda_{js}\) is the wavelength; the ionospheric delays \(I_{r,j}^{s}\) at different frequencies can be expressed as \(I_{r,j}^{s}=\kappa_{j}\cdot I_{r,1}^{s}(\kappa_{j}=\lambda_{j}{}^{2}/\lambda_{1 }{}^{2})\); \(T_{r}^{s}\) represents the slant tropospheric delay; \(e_{r,j}^{s}\) and \(e_{r,j}^{s}\) represent the sum of measurement noise and multipath error for the pseudorange and carrier phase observations. The phase center offsets and variations, tidal loading and phase wind-up are also considered using the existing models with sufficient accuracy ([PERSON] et al., 2017). For the phase center offset (PCO) and phase center variation (PCV) correction of receivers' antennas, the GPS values are used. Because of the different frequencies and signal structure of the individual GNSS, the code biases for the four systems are different from each other in one multi-GNSS receiver. Thus, both inter-system bias (ISB) and inter-frequency bias (IFB) must be taken into account in a combined processing of multi-GNSS observations. The code bias for GPS satellites is set to zero in order to eliminate the singularity between receiver clock and code bias parameters. In other words, all the estimated code biases for the other systems are relative to the biases of the GPS satellites.