## 2 Materials and Methods
### GNSS Tomographic Algorithm  
The current study adopts the tomographic algorithm proposed by [17] (from now on, refer to MM21), with improvements reported in [18; 19]. The advantages of MM21 include its simplicity of implementation, the absence of additional equations for imposing vertical and horizontal constraints, and the fact that there is no need for external observations for initialization. The tomography grid was configured with (7 \(\times\) 10 \(\times\) 20) 1400 voxels, featuring a horizontal resolution of 9 \(\times\) 9 km\({}^{2}\) and a vertical resolution of 500 m (up to 10 km elevation). Figure 1 shows the tomographic grid footprint centered around Houston City. The tomographic domain was rotated by 30\({}^{\circ}\) to include more GPS stations and avoid a large number of empty voxels. The tomographic grid was configured with a vertical terrain-following coordinate [61]. This method transforms the vertical coordinate system to follow the contours of the Earth's surface, ensuring that the vertical levels remain closer to the actual topography. \begin{table}
\begin{tabular}{c c} \hline \hline
**Name of Strategy** & **Strategy Setting** \\ \hline Approach & DD (Double Difference) \\ Choice of observable & LC (Linear Combination of L1 and L2) \\ Cut-off elevation angle & 5\({}^{\circ}\) \\ Sampling interval & 30 s \\ Dry a priori model & SAAS (Saastromainen model) \\ Dry and Wet mapping function & VMF3 \\ Ocean tidal model & FES2004 \\ Solid tide models & IERS10 \\ Orbit and clock & IGS final products \\ Satellite and antenna phase center & IGS14 ANTEX file \\ Tropospheric Gradient Estimation & 2 h \\ Tropospheric delay & Random walk with process noise of 0.02 m\(\sqrt{\text{h}}\) \\ Ionosphere effect & 1-order (LC), 2- and 3-order correction \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Specific processing strategies used in the GAMIT/GLOBK (GG) package. To overcome the matrix rank problem, new slant observations at different azimuths and elevations are computed by least-squares fitting a hypersurface to the SIWV observations [17], given by  
\[\text{SIWV}_{i,k}=a_{i}\text{cos}(\alpha_{i,k})+b_{i}\text{sin}(\alpha_{i,k})+ \frac{c_{i}}{\text{sin}(\epsilon_{i,k})}+d_{i}\ (i=1,\ldots,n;k=1,\ldots,N) \tag{11}\]  
where SIWV\({}_{i,k}\) is the \(k(=1,\ldots,N)\) slant observation at station \(i\) within a specific time window, \(\alpha_{i,k}\) and \(\epsilon_{i,k}\) are the corresponding azimuth and elevation, and \(N\) and \(n\) are the number of slants (per time window) and stations, respectively. The coefficients \(a_{i},\ b_{i},\ c_{i},\ \text{and}\ d_{i}\ (i=1,\ldots,n)\) are then used to create new slants at specific azimuths and elevations, as  
\[\text{SIWV}_{j}=a_{i}\text{cos}\big{(}\alpha_{j}\big{)}+b_{i}\text{sin}\big{(} \alpha_{j}\big{)}+\frac{c_{i}}{\text{sin}\big{(}\epsilon_{j}\big{)}}+d_{i}\ (i=1,\ldots,n;j=1,\ldots,J) \tag{12}\]  
where \(J\) is the number of targets (number of voxels); in our case, targeting the voxel centers from each station location is enough to obtain 100% rank. To exclude the slant contribution of the rays that enter the tomographic domain through its lateral boundaries, the SIWV\({}_{j}\) observation must be corrected. We used the equation proposed by [19] (an improvement of the original equation proposed by [17]), given by  
\[\text{SIWV}^{\prime}_{j}=\text{SIWV}_{j}\text{:}e^{-z/H_{v}} \tag{13}\]  
where SIWV\({}^{\prime}_{j}\) is the corrected SIWV\({}_{j}\), \(z\) is the height at which the rays enter the domain, and \(H_{v}\) is the scale height of water vapor density. This formulation assumes that the background atmosphere is characterized by exponential decaying of the water vapor density. It was found in [19] that keeping \(H_{v}\) constant (normally \(H_{v}=2000\)[62]) does not produce the best tomographic results. Indeed, \(H_{v}\) exhibits both seasonal and diurnal cycles, which influence the tomographic solution.