###### Abstract
[29] pointed out that the integration of GNSS and GRACE produces a higher spatiotemporal resolution than does GRACE alone. The GNSS has the advantages of recovering the TWS with a daily temporal resolution and near-real time, i.e., increasing attention is being paid to high-frequency extreme hydrological events, such as hurricane Harvey along the Gulf Coast [33], and heavy precipitation in Yunnan [17] induced TWS changes. In many studies, the conventional and Laplacian constraints are commonly applied to the GNSS inversion. They both help to solve ill-posed equations when the number of pixels is larger than the number of GNSS stations. In addition, the Laplacian constraints are effective tools used to ensure that the TWS changes vary smoothly between neighboring pixels. The principle of Laplacian constraints is illustrated in Figure 1. There are two kinds of Laplacian constraints. To constrain a target grid point, the simple Laplacian constraints use its four neighboring grid points, the left one, the right one, and the ones on the top and bottom, to force their average to be zero. The complex constraints use all eight neighboring grid points. In essence, the Laplacian constraints require a complex and intricate process to determine the optimal smoothing factor in order to balance the roughness of the TWS in the spatial distribution and the misfit between the model's prediction and observation. Sometimes, determining the optimal smoothing factor using the L-curves [34] or generalized cross-validation (GCV) methods is difficult [35]. Other researchers have attempted to establish different constraint strategies. [PERSON] et al. presented the _a priori_ constraint based on the hydrological models, and the optimal regularization parameter was determined through the iterative least squares estimation instead of the L-curves or GCV method; simulations show that TWS variations from this new constraint can produce more robust inversion results than the traditional Laplacian constraint [20]. [PERSON] et al. proposed a boundary-included inversion algorithm, which reduced the variance by 10%-20% in the boundary regions [12]. However, both the _a priori_ constraint and boundary included inversion methods need to determine the optimal regularization parameters. In addition, the Laplace matrix constitutes tightly fixed constraints, which may result in overly smoothed patterns. To address these issues, we propose a loose spatial constraints strategy. As shown in Figure 1c, the neighboring pixels of a certain pixel are defined as all the pixels within a given radius instead of fixed points. In addition, the TWS at the target point is allowed to change within a threshold value. This threshold value is determined based on a presumed difference between the target point and its surrounding area. These two factors, the given radius and the threshold value, are the key features of the spatial constraints. We aim to improve the spatial distributions of TWS variations across the CONUS in this study. We establish more flexible constraint schemes based on the spatial constraints, which are exempt from the precise determination of the optimal smoothing factor. We first conduct simulations to demonstrate that the inversion with spatial constraints improves the spatial pattern of the TWS. Then, we compare the TWS from different inversions to the TWS provided by the models and time-variable gravity.