## 3 Methodology  
### GNSS Inversion Method  
According to the elastic deformation theory [5; 28; 56; 57], the relationship between vertical displacements and TWS in the EWH reads:  
\[L+v=Gx \tag{4}\]  
where \(L\) represents the vertical displacements measured by GNSS, \(v\) represents the observation residual, \(G\) represents [PERSON]'s function calculated based on the disc model [7; 8], and \(x\) represents TWS to be estimated for each grid cell; the study area is divided into a \(1^{\circ}\times 1^{\circ}\) grid cell. In order to mitigate the edge effects, we extend the study area by 2.5\({}^{\circ}\) and 0.5\({}^{\circ}\) beyond the land and ocean, respectively. The number of grid cells is less than the number of GNSS sites in this study; therefore, the inversion solutions are overdetermined problems. Equation (4) is subject to this condition:  
\[\|Pv\|^{2}\to min \tag{5}\]  
where \(P\) represents the diagonal weight matrix, which is constructed using the inverse of the variance in GNSS vertical displacements. Equation (5) is ill-posed when the number of grid cells is larger than the number of GNSS stations; the Tikhonov regularization method is commonly used for solving ill-posed equations [58]. In addition, to ensure smooth spatial variations in the TWS between adjacent grid cells, a regularization term like Tikhonov is added to Equation (4). The constraint based on the Laplacian operator is as follows [59]:  
\[\lambda Dx=0 \tag{6}\]  
where \(\lambda\) represents the smoothing factor that adjusts the relative weight between the model misfit and roughness. \(D\) represents the Laplacian operator with the following kernels: \(\left[ \begin{array}{ccc}0&1&0\\ 1&-4&1\\ 0&1&0\end{array} \right]\). Thus, Equation (5) is written as follows after adding the Laplacian constraints:  
\[\|Pv\|^{2}+\lambda^{2}\|Dx\|^{2}\to min \tag{7}\]  
where the first term \(\|Pv\|^{2}\) represents the model misfit, and the second term \(\|Dx\|^{2}\) represents the roughness of the solution. The TWS is estimated as follows:  
\[x=\ (G^{T}PG+\lambda D^{T}D)^{-1}G^{T}PL \tag{8}\]  
In addition to the Laplacian constraint strategy mentioned above, the conventional constraint strategy is also applied in this study, which is similar to the Laplacian constraints in Equation (6), but the Laplacian operator is replaced with the identity matrix. The conventional constraint strategy addresses the issue of ill-posed equations but does not consider the smoothness of the TWS changes between adjacent grid cells. The formula is as follows:  
\[\lambda Ix=0 \tag{9}\]  
where \(I\) represents the identity matrix. The TWS is inverted as follows:  
\[x=\ (G^{T}PG+\lambda I)^{-1}G^{T}PL \tag{10}\]The optimal smoothing factors need to be determined before inversion for the conventional constraints and Laplacian constraints. We use the GCV approach to estimate the optimal parameters. We utilize the GNSS vertical displacements from 1 January 2018 to 31 August 2022 at 1437 GNSS stations across the CONUS and run the inversion with various values of \(\lambda\). Subsequently, we obtain the curves between the smoothing factors and the sum of squared residuals obtained through cross-validation (CVSS). As shown in Figure 6, the optimal parameters are 3.2 for conventional constraints and 1.3 for Laplacian constraints, as in these cases, the CVSS values are minimized. As mentioned above, the principle of spatial constraints is as follows: the absolute difference between the EWH of the \(k\)th grid cell and the average EWH of all the grid points in the circle with a radius of \(r\) centered on the \(k\)th grid cell is within a certain range. This is illustrated by the following inequality:  
\[|x_{k}-\frac{1}{n}(x_{k_{1}}+x_{k_{2}}+\cdots+x_{k_{n}})|\leq e \tag{11}\]  
where \(x_{k}\) represents the EWH of the \(k\)th grid cell; \(x_{k_{1}},\ x_{k_{2}},\ \cdots,\ x_{k_{n}}\) represent the EWH of all the grid points within a circle of radius \(r\) centered on the \(k\)th grid cell; and \(e\) represents the set threshold value.