## 3 Methodology
### Closed-Loop Simulation  
We conduct a closed-loop simulation to assess the performance of spatial constraints for recovering the TWS variations. The procedure is described as follows: Firstly, we extract soil moisture and snow water data from the monthly GLDAS \(1^{\circ}\times 1^{\circ}\) dataset in 2020 across the CONUS, which is used as the input signal. Next, we forward calculate the monthly vertical displacements at each GNSS station. The noise signals at each GNSS station are obtained by averaging the STD of vertical displacements from 1 January 2020 to 31 December 2020 in each month. Then, we apply different constraint strategies to invert the monthly spatial distributions of the TWS and obtain the TWS time series. Finally, we calculate the amplitudes and STDs of monthly TWS from the conventional (CC-TWS), Laplacian (LC-TWS), and spatial constraints (SC-TWS).  
For the spatial constraint inversion algorithm, two empirical parameters, the radius \(r\) and the threshold \(e\), jointly control the spatial smoothness of the TWS variations, which, in turn, affects the oscillation amplitude of the TWS time series. We illustrate the performance of the spatial constraints in TWS inversion using different combinations of \(r\) and \(e\) (see Table 2). The amplitudes and average STDs are shown in Table 2, the average STDs of the SC-TWS are lower than those of the CC-TWS and LC-TWS, indicating that the spatial constraint inversion results are closer to the input TWS. The conventional constraints exhibit the minimum oscillation amplitude, and the spatial constraints show the maximum oscillation amplitude.  
We use the following equation to quantify the improvement of TWS using the spatial constraint method:  
\[\Delta\text{STD}=\frac{\text{STD}_{\text{CC/LC}}-\text{STD}_{\text{SC}}}{ \text{STD}_{\text{CC/LC}}}\times 100\% \tag{16}\]  
where \(\text{STD}_{\text{CC/LC}}\) represents the average STDs of the conventional or Laplacian constraints, \(\text{STD}_{\text{SC}}\) represents the average STDs of the spatial constraints, and \(\Delta\)STD represents the improvement of the spatial constraint method. From Table 2, we can see that the spa  
\begin{table}
\begin{tabular}{c c c} \hline \hline
**Constraint Strategy** & **Amplitude (mm)** & **Average STDs (mm)** \\ \hline input & 52.59 & \\ CC & 50.06 & 14 \\ LC & 50.21 & 14 \\ SC15045 & 51.03 & 11 \\ SC20045 & 50.76 & 12 \\ SC25055 & 50.74 & 11 \\ SC30060 & 50.82 & 11 \\ SC35065 & 50.76 & 11 \\ SC40060 & 50.68 & 13 \\ SC45075 & 51.00 & 12 \\ SC50075 & 50.89 & 12 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Annual amplitudes of input GLDAS and recovered TWS and average STDs between input GLDAS and recovered TWS across the CONUS in 2020. CC and LC represent conventional and Laplacian constraints, respectively. The first three digits and the last two digits after the letter SC represent the values of the parameters \(r\) and \(e\) of the spatial constraint inversion method, respectively. For example, SC15045 represents the parameters \(r=150\) km, \(e=45\) mm, and so on.  
tial constraint inversion is improved by 7-21% compared with the inversion using the conventional and Laplacian constraints.