###### Abstract  
Global Navigation Satellite System (GNSS) vertical displacements are widely used to infer terrestrial water storage (TWS) variations. The traditional Laplacian inversion requires dedicated efforts to determine the optimal parameters, which has an important effect on the spatial patterns. In this study, we develop a new GNSS inversion method with flexible spatial constraints. One major merit is that the new method only requires loose boundary conditions rather than optimal parameters. A closed-loop simulation shows that the inversion using spatial constraints is improved by 7-21% compared with the Laplacian constraints. We apply this method to 18 watersheds across the Contiguous United States (CONUS) to infer daily TWS variations from January 2018 to August 2022. The results show that the amplitudes of monthly TWS time series from the spatial and Laplacian constraints are comparable to the Gravity Recovery and Climate Experiment (GRACE) Follow-On (GFO) in 16 watersheds. Furthermore, the standard deviation between the spatial constraints and GFO is at the same level as that between the Laplacian constraints and GFO. We also extract the daily TWS variations caused by heavy precipitation events in California. Our results demonstrate that spatial constraint inversion supplements the existing constraint strategies of GNSS inversion in hydrogeodesy; therefore, spatial constraint inversion can be an alternative tool for GNSS inversion. t 202315753 10.3390/rs15245753 10.3390/10  
period in the western United States, which shows consistency with the decreases in precipitation and streamflow [6]. [PERSON] et al. recovered the seasonal TWS variations in the northwest of the CONUS; they show that the GNSS can be independently used to estimate the TWS variations. The mass changes in dry and wet periods can be inferred using the GNSS, and it can fill the gap between GRACE and GFO [1]. [PERSON] and [PERSON] estimated the monthly evapotranspiration variation using the GNSS-derived TWS in the Mississippi River basin; the results indicate strong agreement between the evapotranspiration variations obtained from the GNSS and other remote sensing datasets [9]. [PERSON] et al. recovered the TWS in Southwest China using an independent component analysis (ICA)-based inversion method; the results show that the TWS variations from ICA-based inversion highly agree with the GRACE and hydrologic models, especially in the northwest of the study region with an extreme and complex terrain, where the TWS derived from the original GNSS data was overestimated [19]. In areas with a sparse distribution of GNSS stations, the Slepian basis function (SBF) is another effective tool used for TWS inversion. [PERSON] and [PERSON] proposed the SBF method for surface mass variation inversion in the Australian continent for the first time using GNSS vertical displacements [21]. Subsequently, [PERSON] et al. inferred the TWS variations in mainland China and Brazil; they developed a novel drought index and assessed the severity of hydrological drought [22; 23]. [PERSON] et al. developed a hydrological drought index at multiple scales based on the TWS using the SBF and found that it is more robust in quantifying Brazil's hydrological drought [24]. [PERSON] et al. inverted the TWS variations in southwest China; they found that the amplitudes of the TWS using the SBF are larger than those of GRACE/GFO, and the filter radius of the SBF inversion method can be determined according to the average distance between the GNSS stations [25]. The TWS inferred using the Green's function and SBF methods show good agreement with the GRACE/GFO and hydrologic models, but there are differences in the spatial and temporal patterns in some regions [26]. Other recent studies focused on obtaining more reliable TWS estimates through the joint inversion of the GNSS and GRACE/GFO in the CONUS [27; 28; 29] and China [26; 30; 31; 32]. [PERSON] et al. [27] and [PERSON] et al.