AbstractGridded terrestrial water storage (TWS) variations observed by GRACE or GRACE-FO typically show a spatial correlation structure that is both anisotropic (direction-dependent) and non-homogeneous (latitude-dependent). We introduce a new correlation model to represent this structure. This correlation model allows GRACE and GRACE-FO data users to get realistic correlations of the TWS grids without the need to derive them from the formal spherical harmonic uncertainties. Further, we found that the modelled correlations fit the spatial structure of uncertainties to a greater extent in a simulation environment. The model is based on a direction-dependent Bessel function of the first kind which allows to model the longer correlation lengths in the longitudinal direction via a shape parameter, and also to account for residual GRACE striping errors that might remain after spatial filtering. The global scale and shape parameters vary with latitude by means of even Legendre polynomials. The correlation between two points transformed to covariance by scaling with the standard deviations of each point. The covariance model is valid on the sphere which is empirically verified with a Monte-Carlo approach. The covariance model is subsequently applied to 5 years of simulated GRACE-FO data which allow for immediate validation with true uncertainties from the differences between the input mass signal and the recovered gravity fields. Four different realisations of the point standard deviations were tested: two based on the formal errors provided with the simulated Stokes coefficients, and two based on empirical standard deviations, where the first is spatially variant and temporally invariant, and the second spatially invariant and temporally variant. These four different covariance models are applied to compute TWS time series uncertainties for both the fifty largest discharge basins and regular grid cells over the continents. These four models are compared with the true uncertainties available in the simulations. The two empirically-based covariance models provide more realistic TWS uncertainties than the ones based on the formal errors. Especially, the empirically-based covariance models are better in reflecting the spatial pattern of the uncertainties of the simulated GRACE-FO data including their latitude dependence. However, these modelled uncertainties are in general too large. But with only one global scaling factor, a statistical test confirms the equivalence between the empirically-based covariance model with temporally variable point standard deviations and the true uncertainties. Thus at the end, this covariance model represents the closest fit in the simulation environment. The simulated GRACE-FO data are assumed to be very realistic which is why we recommend the new covariance model to be further investigated for the characterisation of real GRACE and GRACE-FO terrestrial water storage data.