Evaluating long-term water storage trends in small catchments and aquifers from a joint inversion of 20 years of GRA CE/GRA CE-FO mission data

More than 20 yr of measurement data of the gravity missions GRACE (Gravity Recovery And Climate Experiment) and GRA CE-FO (GRA CE-Follo w-On) allo w detailed in vestiga-tions of long-term trends in continental terrestrial water storage (TWS). Ho wever , the spatial resolution of conventional GRACE/GRACE-FO data products is limited to a few hundred kilometres which restrains from investigating hydrological trends at smaller spatial scales. In this study GRACE and GRACE-FO data have been used to calculate TWS trends with maximized spatial resolution. Con ventionally, GRA CE/GRA CE-FO is presented as a series of either unconstrained gravity ﬁelds post-processed with spatial low pass ﬁlters or constrained inversions commonly known as Mascon products. This paper demonstrates that both approaches to suppress spatially correlated noise are mathematically equivalent. Moreover, we demonstrate that readily inverting all available sensor data from GRA CE/GRA CE-FO for a single TWS trend map, together with annual variations and a mean gravity ﬁeld, provides additional spatial detail not accessible from the standard products. The variable trade-off between spatial and temporal resolution as a unique feature of satellite g ravimetr y allows for gravity products that are tailored towards speciﬁc geophysical applications. We show additional signal content in ter ms of long-ter m water storage trends for four dedicated examples (Lake Victoria, North-west India, Bugachany Reservoir and High Plains Aquifer) for which external information from other remote sensing instruments corroborates the enhanced spatial resolution of the new mean-ﬁeld trend product.


I N T RO D U C T I O N
Dedicated satellite missions that enable the precise and continuous surv e y of spacecraft trajectories in low altitude orbits around the Earth have revolutionized the ability to map the Earth's global gravity field and its temporal variations (Flechtner et al. 2021 ).The Gravity Recovery and Climate Experiment (GRACE, Tapley et al. 2019 ) and its successor GRA CE Follo w-On (GRA CE-FO , Landerer et al. 2020 ) allow the calculation of a new gravity field model from sensor data accumulated over only 30 d.The resulting time-series e xtends ov er more than 20 yr and has spurred entirel y ne w and innov ati ve applications of satellite g ravimetr y in many branches of the physical Earth sciences.GRACE and GRACE-FO data allowed for the first time the unambiguous mapping of large-scale mass losses from continental ice sheets (Shepherd et al. 2012 ) and mountain glaciers (Cirac ì et al. 2020 ).The missions also provide direct estimates of the increase in ocean mass due to enhanced meltwater inflow as a consequence of ice mass loss (Cazenave & Moreira 2020 ), which allows for the separation of sea-level rise solel y caused b y ocean w arming as a direct measure of the Earth's energy imbalance (von Schuckmann et al. 2016 ).Satellite gravimetry has also triggered new research in the solid Earth sciences by providing new observational insight into glacial isostatic adjustment (Caron & Ivins 2020 ), the temporal evolution of megathrust earthquakes (Panet et al. 2022 ), and ongoing tectonic processes in active continental collision zones like the Tibetan Plateau (Braitenberg & Shum 2017 ).
In addition, the data from GRACE and GRACE-FO have significantly impacted h ydrology, h ydrometeorology and the climate sciences.For the first time, variations in terrestrial water storage 1002 C The Author(s) 2023.Published by Oxford University Press on behalf of The Royal Astronomical Society.This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( https://creati vecommons.org/licenses/by/4.0/ ), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.Downloaded from https://academic.oup.com/gji/article/236/2/1002/7459368 by GFZ Potsdam user on 09 January 2024 Evaluating long-term water storage trends 1003 (TWS) in all possible storages from the surface down to the deepest aquifers are mapped over all continental areas of the world.Besides identifying seasonal variations in TWS (Tapley et al. 2004 ), the satellite record provides insight into interannual variations related to large-scale climate modes (Pfeffer et al. 2022 ), the characteristics of hydrometeorological extremes like extended periods of drought (Thomas et al. 2014 ) or heavy flooding events (Chen et al. 2010 ), as well as the quantification of atmospheric net-water fluxes consisting of precipitation and e v apotranspiration (Eicker et al. 2020 ).GRA CE data also allo wed for the first time to identify large-scale groundwater losses (Famiglietti 2014 ) due to extensi ve anthropo genic w ater extraction for irrigation.Such human interventions into the terrestrial water cycle are expected to increase with increasing surface temperatures due to continuously rising greenhouse gas concentrations in the atmosphere.Even today, contemporary trends emerging in the GRACE record align well with long-term projections from global coupled climate models (Jensen et al. 2019 ), thereby nicely corroborating inferences obtained from such models for the future evolution of our planet.The importance of continued monitoring of TWS from space utilizing satellites has also been acknowledged by the Global Climate Observing System of the World Meteorological Organization by listing TWS as a new Essential Climate Variable (ECV) in its 2022 Implementation Plan (World Meteorological Organization 2022 ).
In contrast to the gravity gradiometer mission GOCE that directly observes functionals of the Earth's gravity field along the satellite orbit (Rummel et al. 2011 ), the spacecrafts of each of the twin-satellite missions GRACE and GRACE-FO instead act as proof masses in Earth's gravitational field.Precise monitoring of the relative distance, velocity and accelerations between the two identical satellites trailing each other approximately 200 km apart allows for inferring the gravitational disturbances of the mass distribution at or underneath the Earth's surface.The tracking is realized with microwave frequencies in the K-and Ka-bands for both missions, and on GRACE-FO additionally with an innov ati ve laser interferometr y instr ument (Abich et al. 2019 ) that of fers e ven higher tracking precision.Nevertheless, satellite measurements require a mathematical inversion to obtain surface mass changes at Earth's surface.
This inversion poses several challenges which directly impact the quality of terrestrial water storage estimates derived from satellite g ravimetr y.Each satellite is affected by the gravitational acceleration changes at satellite altitude, which is smoothed compared to the gravitational attraction at Earth's surface due to the large distance of the satellite to the field-inducing masses.We note that lo wer orbits w ould increase orbit perturbations and thus the sensitivity of gravity missions.Ho wever , higher atmospheric densities and, consequently, more considerable atmospheric drag acting on the satellites would require additional propulsion energy to compensate for the loss of momentum and subsequently orbit altitude.Therefore, the initial orbit height of 500 km is a carefully chosen compromise between mission duration and mission sensitivity.
During the inversion, the observed signal needs to be propagated from satellite altitude to Earth's surface, to determine water storage changes.This downward continuation process also amplifies measurement noise depending on the spatial scale.Smaller scales are subject to stronger amplification, resulting in well-known high-frequency spatial noise patterns (e.g.Ilk et al. 2002 ).Two approaches for noise reduction can be distinguished: Traditionally, gravity field parameters are solved in terms of (unconstrained) spherical harmonics that are subsequently smoothed in space to diminish the impact of noise contained in the poorly determined smaller spatial scales (Wahr et al. 1998 ).Spatial smoothing might also be performed with anisotropic filters objecti vel y deri ved from the characteristics of the error covariance matrix of a typical monthly solution of the global gravity field derived from GRACE data (Kusche 2007 ).The second approach includes a priori constraints about surface mass variability for geographic blocks (or spherical caps) to stabilize the inversion directly without requiring subsequent post-processing.Pioneered for GRA CE/GRA CE-FO applications by Luthcke et al. ( 2013 ), Mascon solutions are now available from many groups that also include the GRA CE/GRA CE-FO Science Data System centres Center for Space Research at University of Austin (CSR, Save et al. 2016 ) and Jet Propulsion Laboratory (JPL, Watkins et al. 2015 ) as well as from the Goddard Space Flight Center (GSFC, Loomis et al. 2019 ).A comprehensiv e ov erview on the development of Mascon solutions can be found in Antoni ( 2022 ).
Both approaches have the same effect by acting as a spatial lowpass filter to reduce the noise introduced by the inversion process and are mathematically equivalent.That is evident mainly because the constrained inversion can be rewritten as a filter matrix acting on the unconstrained solution.For further consideration, we treat both approaches as a low-pass filter acting on a noisy input surface mass change field.A key design parameter in both approaches is the signal-to-noise ratio (SNR) of the final mass change fields.In other words, one has to decide how strong the implicit or explicit low-pass filter should be to obtain a data product with a suitable noise level.This decision can be based on external models for the observation noise and expected signal.Regardless of how the filter is derived, it is beneficial to filter as little as possible since additional smoothing will reduce the ef fecti ve spatial resolution of the resulting mass grids as not only noise but also signal is filtered (e.g.Vishwakarma et al. 2017 ;Wiese et al. 2016 ).Additionall y, it will ine vitabl y enhance the adverse ef fects of spatial leakage, which is one of the most considerable obstacles in applying GRACE data for scientists with a non-geodetic background.
The overall noise level in the unconstrained surface mass change estimates is a significant factor in determining the filter strength.By accumulating observations of a more extended period than 30 d, the random noise components tend to average out, thus reducing the noise level in the derived unconstrained gravity field and the subsequently inverted surface mass estimates.In other words, we can directly benefit from the reduced noise le vel b y appl ying a weaker filter, thus retaining more signals with less spatial leakage so that processes at smaller spatial scales might be studied.Ef fecti vel y, spatial resolution is traded against the temporal resolution.
In this contribution, we contrast long-term water storage trends estimated from monthly GRA CE and GRA CE-FO solutions against trends based on accumulated observations over the same period.In Section 2 , we show the equi v alence of post-processing filters and constrained inversion, detail the deri v ation of our long-term TWS trend estimate and outline the data sets used in this study.Section 3 contains an intercomparison of trend signals derived from different data sets in small regions, as well as a comparison with independent observ ations.Finall y, we discuss the advantages and disadvantages of long-term water storage trend estimates and highlight areas of use where this complementary data set offers benefits over the more traditional analysis of monthly GRA CE/GRA CE-FO water storage products.

D ATA A N D M E T H O D S 2.1 A conceptual look at water storage trend estimates based on GRA CE/GRA CE-FO data
For better accessibility, TWS anomalies based on GRACE and GRACE-FO data are typically given in equivalent water height (EWH) on geog raphic g rids.As obser v ation-deri ved data sets, the treatment of the (inevitable) measurement noise has a significant impact on the quality of the TWS data products, which are either based on the output of a constrained inversion (CI), as in the popular Mascon products, or derived from post-processing unconstrained gravity field solutions b y appl ying spatial low-pass filters (SLPF).In practice, different strategies exist for the determination of the filter parameters.The DDK filter of Kusche et al. ( 2009 ) uses a synthetic noise model for GRACE observations based on a single month's satellite ground track coverage.That works very well for most of the GRA CE/GRA CE-FO period, except for months where the satellites exhibit a repeat orbit which deteriorates the filtered surface mass change fields.Horvath et al. ( 2018 ) designed the VDK filter to additionally consider changes in the orbit geometry and the actual measurement noise le vel, which e ventuall y increases during a mission due to energy limitations leading to reduced thermal stability, and other hardware degradations.Filtering with VDK is, therefore, directly comparable to a constraint inversion based on the same data.The signal covariance model of both the DDK and VDK filters is based on an isotropic po wer la w function.That means the filter expects the same signal magnitude and spatial correlation at each point on Earth's surface.Compared to the VDK filter, Mascon solutions onl y dif fer in complexity and deri v ation of the signal covariance model, which in their case may be variable in space and time.A special case are the JPL Mascons, which also consider temporal correlations between adjacent months.An advantage of the Mascons is their parametrization in the space domain, which makes formulating complex spatially varying constraints easier than the otherwise utilized spherical har monics.Never theless, it is mathematically possible to introduce similar information in the spectral domain (Chen et al. 2021 ).
Even though the CI and SLPF approaches seem very different, they are closely related.Both approaches aim at suppressing spatially correlated noise, which renders the usage of unconstrained or unfiltered GRA CE/GRA CE-FO EWH grids impossible.This noise increases inversely with the spatial scale, so CI and SLPF approaches act as low-pass filters.The main difference between the approaches is where the low-pass filter effect comes into play.In CI, constraints in the form of a signal variance matrix are applied directly in the inversion process to mitigate the effect of measurement noise on the water storage estimates.Conversely, in SLPF, a spatial low-pass filter is applied to unconstrained GRACE/GRACE-FO gravity field solutions before propagating to EWH grids.The final results of both approaches are comparable and mathematically interchangeable.
As mentioned abov e, ev ery constraint can be expressed as a filtering operation, which is also the idea behind the widely used DDK filter (Kusche et al. 2009 ).To show the equi v alence between constraint inversion and applying a spatial filter, we start with the least squares adjustment (LSA) to estimate the unconstrained/unfiltered solution ˆ x U .Given the satellite measurements l , the linearized functional model A , and the residuals e , the observation equations for this LSA problem are given by l = Ax U + e . (1) The solution to this system of equations is obtained by forming the normal equation matrix N and right-hand side n , where P is the inverse v ariance-cov ariance matrix of the observation noise.Constraints can be applied by adding the inverse of the signal covariance matrix S .Thus, the constrained solution can be written as where N ˆ x U = n is the system of normal equations of the unconstrained estimate.The spatial filter matrix equi v alent to the constrained estimate is then defined as which can be rewritten to using the Woodbury matrix identity (Woodbury 1950 ).In eq. ( 5), ˆ x = N + is the noise covariance matrix of the unconstrained solution and S is the expected signal covariance matrix.The filter or the constrained adjustment is, therefore, fully described by the signal and noise covariance models S and ˆ x .In the case of monthly total water storage products, the signal covariance matrix S should reflect the month-to-month variability of water storage and ˆ x is a measure of the observation noise present in an unconstrained monthly solution.Note that these considerations are independent of the chosen gravity field parametrization, so they apply to both spherical harmonics-and Mascon-based solutions.
We can use eq.( 5) to investigate how the filter properties change with different noise levels.For low noise ( ˆ x → 0 ), the filter matrix tends towards the identity matrix, resulting in a weak filter.That is also the case if the expected signal S is drastically larger than the noise.For high noise levels, W will tend towards zero, which means both the signal and the noise are aggressi vel y filtered.
Ideally, a post-processing filter maximizes the SNR in the soughtafter water storage grids by removing as much noise as possible while retaining all the signals.Since applying a filter to the EWH grids will affect both signal and noise, we also want to choose the least aggressive filter, which allows us to reach a suitable SNR.A reasonable goal is, therefore, to look for ways to reduce the input noise, which then allows applying weaker filters and constraints.So, we can reduce the GRA CE/GRA CE-FO measurement noise by sacrificing temporal resolution.That comes somewhat intuitively as restricting the possible temporal variation acts as an averaging process in time.Applying this to our previous considerations means that reducing the spatial resolution and, thus, reducing the measurement noise, the filters we need to apply can be less aggressive.
The conclusion we can draw here is that filters and constraints tailored to month-to-month water storage variations have to deal with a higher inherent measurement noise level than filters and constraints for a multiyear water storage trend due to the different averaging periods.Subsequently, trends computed from monthly GRA CE/GRA CE-FO data products are filtered too aggressi vel y.That means crucial trend information is lost, especially in smaller catchments and aquifers.
If we directly formulate constraints which take into account the expected trend signal and the reduced noise level resulting from the long accumulation period, we can arrive at a better water storage trend estimate, which preserves more small-scale features.This approach is followed in long-term trends co-estimated with mean gravity field models like GOCO06s, ITSG-Grace2018s or Tongji-Grace2018k.It is, therefore, reasonable to compare these co-estimated trends with trends based on monthly water storage products.In the following, we want to determine the impact of tailored constraints on the derived water storage trends, compared to qualitative and quantitative trends from monthly data products.

Intercomparison of water storage trends based on monthly data sets
First, we compare different monthly products derived with different approaches.To that end, we estimate long-term trends, together with constant and seasonal terms, from the RL06 Mascons solutions of the Center of Space Research (CSR, Center for Space Research 2020 ), the Jet Propulsion Laboratory (JPL, NASA/JPL 2019 ), and the Goddard Space Flight Center (GSFC, Loomis et al. 2019 ), as well as post-processed spherical-harmonics-based solutions of the Inter national Combination Ser vice for T ime-variable Gra vity Fields (COST-G) and the German Research Centre for Geosciences (GFZ).The latter two are provided via the GFZ-hosted GravIS portal (Gravity Infor mation Ser vice) 1 as g ridded water storage anomalies.All these data sets have the same signal definition on the continents and are, therefore, directly comparab le w hen looking at water storage variations.Fig. 1 shows the estimated trends from 2002-04 to 2021-01 for each of the products in the comparison.We can see the most prominent differences between spherical-harmonics-based and Mascon-based solutions in regions where large horizontal gradients in water storage trend values occur.This mainly concerns regions with significant ice mass loss, such as Greenland, Antarctica, the Gulf of Alaska and other coastal areas.Within the continental areas, all the products are similar, with areas in Nor ther n India and China as notable exceptions.
Comparing Mascon solutions on a pix el-per-pix el basis without interpolation is difficult because of the different Mascon placements and shapes.We can, ho wever , evaluate the spherical-harmonicsbased solutions on an arbitrar y g rid and average the water storage values in each Mascon cell.That allows us to directly compare spherical harmonic and Mascon solutions and perform indirect inter-comparisons between the different Mascon products.Fig. 2 shows the differences between the two spherical-harmonic-based solutions used in this study, and each of the three Mascon solutions.In addition to the regions where we were already able to identify differences in the absolute trend values, we can now also see other areas with minor differences.The water storage trends of the CSR and JPL Mascons primarily differ in spatially confined regions like Nor ther n India, China and around Lake Victoria.That can be attributed to the inhomogeneous signal covariance model used in the Mascon process compared to the isotropic one used in the VDK filter.The GSFC Mascons fit closest to the spherical-harmonics-based solutions with differences below 10 mm yr −1 EWH on the continents.Fur ther, we can obser ve that the two spherical-har monics-based solutions are very similar except for slightly reduced striping artefacts in the COST-G products.Overall, all products show a close agreement on the continents.

Trend co-estimated with mean gravity field
To calculate mean field trend (MFT) trend estimate, we take the system of normal equations from ITSG-Grace2018s (May er -G ürr et al. 2018 ; Kvas et al. 2019 ), which consists of a static constituent, trend 1 http://gravis.gfz-potsdam.de/home, visited 2023-01-04 and annual oscillation.We combine it with normal equations assemb led from appro ximately two and a half years for GRACE-FO GPS tracking and microwave ranging data.The parametrization of the individual constituents is the same as for the GRACE contribution, which models the static component of Earth's gravitational field in spherical harmonics up to degree and order 200, and all temporal constituents up to degree and order 120.The two systems of normal equations are then added, a constraint is applied to the annual oscillation, and both static and annual oscillation are subsequently eliminated.That leaves us with an unconstrained system of trend normal equations up to degree and order 120, which implicitly contains static and annual parameters.Before solving the system of normal equations, we transform it from the spherical harmonics domain into the space domain.The reason for this is that regionall y v arying constraints can be easier applied there, because in the spectral domain, they will result in a dense regularization matrix (Kvas et al. 2021 ;Chen et al. 2021 ).Our target parameter space are radial basis functions as defined in Eicker et al. ( 2014 ) with a diagonal regularization matrix.Instead of pre-defining the signal covariance matrix used to regularize the trend, we co-estimate it with the trend parameters using variance component estimation (VCE).That means we have a purely data-driven constraint for the accumulated trend estimate without any prior information on the signal level.
On the other hand, deriving a trend from a time-series of constrained or post-processed solutions is typically done by performing regression in the form of where in addition to the trend parameters, a constant value and an annual oscillation are estimated.In some cases, additional sinusoidal parameters, for example, a 161-d or semi-annual oscillations are set up.This deterministic model can be estimated in a multi v ariate linear regression model of the form where A is the design matrix containing the temporal basis functions, L is the observation matrix which contains the filtered or constrained solutions ˆ x T C, k in each row.The solution to this regression problem is subsequently given by We added the subscripts to the filter matrices because of their time dependency.In the case of GRA CE and GRA CE-FO , the time dependence is a result of varying accuracy due to changing orbit geometry, satellite altitude or instrument quality.Furthermore, some processing centres use a time-varying signal model.For further consideration, we assume an idealized case where we have stationary signal and noise properties.That allows us to find a direct relation between constrained and unconstrained model estimates given by ˆ where W describes a filter with properties representative for the whole time-series.More concisely, we can relate the constrained and unconstrained trend with ˆ a CT ≈ W ˆ a UT .Even though this is a very idealized look at constrained trend estimates, we can conclude that fitting a trend to constrained or post-processed monthly products are roughly equivalent to estimating the trend from unconstrained or unfiltered solution and then applying a filter with the same properties.As a consequence, the resulting trend will have roughly the same properties in terms of signal attenuation and spatial leakage as each input solution.If we revisit eq. ( 5), we can see that the constraints tailored to the monthly signal and noise content represented by W might be a sub-optimal fit for the trend estimate.In regions with large seasonal water storage, the signal covariance S will overestimate the trend.More importantly, the noise level in the unconstrained trend will be significantly lower than in a single monthly solution if the accumulation period is long enough.Thus W will generally result in a too aggressive low-pass filter when long-term trends are concerned.
To make our MFT estimate comparable with the other data sets, we added degree 1 rates obtained from GFZ GravIS2 and removed the effect of GIA using the ICE6G-D model (Peltier et al. 2018 ).The resulting trend propagated to the space domain can be found in Fig. 3 .We can see that most features are sharper and better localized than in the trends based on monthl y solutions.Howe ver, we can also see ringing effects in regions with large signal gradients, like Greenland, Alaska and Antarctica.This is a direct consequence of the maximum spherical har monic deg ree of 120 to which the initial ITSG-Grace2018 and GRACE-FO systems of trend normal equations are developed.While suitable for applications on monthly time scales where smaller spatial scales are smoothed more aggressi vel y due to lower SNR, this chosen maximum degree is not sufficient for long-term trend estimate.Combined with our data-driven regularization technique, which cannot discriminate between geophysical signal and ringing effects, this diminishes the usability of the MFT estimate in regions with larger spatial gradients.These artefacts can be reduced with a higher spherical harmonic maximum degree of the base normal equations, where sharper gradients can be represented or (additionally) a tailored regularization which suppresses these spurious signals.Before the MFT estimate can be released as a robust data product, these points need to be addressed, however, for the study at hand they play a minor, as our goal is to offer a proof of concept that the MFT provides a theoretically higher spatial resolution.

E X A M P L E S O F H Y D RO L O G I C A L T R E N D S I G N A L S
This section will take a closer look at four smaller-scale hydrological trends.These examples will showcase the value of the higher spatial resolution of the MFT compared to trends estimated from monthly solutions.As shown above, the differences between the trends of monthly solutions are relatively small.Thus, we only use the COST-G and the GSFC Mascons trend for comparison.

Lake Victoria
Lake Victoria in Eastern Africa is the second-largest freshwater lake in the world.The lake's w ater le vel has risen since around 2006 about three meters, which induces a strong trend signal in the gravity field.Fig. 4 shows the spatial patterns of the TWS trends around Lake Victoria.The trend is far better localized in the lake for the MFT product than for monthly solutions.For independent comparison, we use a volume change time-series derived from satellite altimetr y obser v ations from DAHITI ( dahiti.dgfi.tum.de, Schw atke et al. 2015 ) and a constant lake surface extent derived from Global Surface Water ( https://global-surface-water.appspot.com/, Pekel et al. 2016 ).The latter is an acceptable simplification since the seasonal and internannual variations of the surface area of Lake Victoria are small compared to its total extent.The altimetric volume change amounts to a total mass difference of 129.6 Gt between April 2002 and January 2021.In order to get the accumulated mass change from the TWS trends, the spatial integration area has to be defined.We chose to integrate first over the lake area (58 000 km 2 ) and then also add a 50 km buffer to the lake shoreline to capture leakage signals (139 000 km 2 ).The buffer width is chosen subjecti vel y, as our primary interest is to gauge the spatial resolution and spatial localization of the different trend products rather than a TWS analysis.For this reason, we also omit uncertainty estimates for the derived TWS trend values.Integrating only over lake area leads to an  underestimation of the mass change in all three trends, with 76.1 Gt for the MFT and 34.9 Gt and 38.8 Gt for the COST-G and GSFC Masons trends, respecti vel y.Taking the 50 km buf fer into account improves all estimates, with 131.3 Gt for the MFT, and 75.2 and 82.3 Gt for COST-G and GSFC Masons trends.Within the extended integration area, the MFT almost captures the altimetry-derived mass change, while the trends derived from monthly solution still exhibit substantial underestimation.That proves the advantage of the significantly better spatial resolution of the MFT in localizing smaller scale hydrological signals.Note that extending the integration area to counteract signal attenuation only works for cases where there is little to no signal outside of the region of interest.Extending the integration region in areas with strong neighbouring signals will mix signal sources and thus produce results that are hard to interpret.Ho wever , in such cases, the MFT still proves advantageous because due to the better spatial localization, as this example showcases, signal cross-contamination will be lo wer , and the magnitude will be higher even if only the catchment boundaries are considered.

Northwestern India
The following example takes a look at the well-known trend signal in Nor thwester n India (e.g.Rodell et al. 2018 ).While the two trends calculated with the monthly data only show a strong trend in the Nor thwester n par t of India, the TWS trend of the MFT solution shows a considerably more detailed spatial pattern than the trends derived from the monthly solutions Fig. 5 .The marked ne gativ e MFT trends to the northwest of the city of Delhi in the Ganges-Indus plains and on the eastern side of the Arav all range largel y coincide with regions that are characterized by irrigation that is to over 85 per cent fed from groundwater, according to the Food and Agricultural Organization (FAO, https: //data.apps.fao.org/aquamaps/, see outlines in Fig. 5 ).The ne gativ e TWS trends can thus be related to groundwater overuse as stated in previous studies, but are spatially attributed in a much better way to the irrigation areas in the MFT data set.Another smaller area with ne gativ e TWS trends in MFT in the Upper Indus Valley is most probably caused by melting high mountain glaciers in this region (Brun et al. 2017 ).All these spatially detailed trends cannot be distinguished in the trend maps of the two monthly products.

Boguchany reservoir
The Bo guchan y Reservoir dams the Angara River in Siberia near the town of Kodinsk.Though the works on the dam started as early as 1974, the construction w as onl y finished in 2012, after which the reservoir was filled until 2016.The official capacity of the reservoir amounts to 58.2 km 3 (SE Solutions & Ecoline Environmental Assessment Centre 2006 ) and it extends from the dam in Kodinsk to the Ust-Ilimsk dam located 357 km upstream with a width ranging from 1.2 to 15 km (Jagus & Rzetala 2013 ).Fig. 6 shows the TWS trend signals around the reservoir, again the MFT and the two trends derived from the monthly solutions of COST-G and GSFC Mascons.While the MFT shows a clear and robust positive signal in the reservoir's area, the two monthl y-solutions-deri ved trends display only a weak signal.We first compute the mass change for an area that represents the reservoir outlines by adding a 15 km buffer to the corresponding river segment.Here, all trend products underestimate the reservoir capacity with 15.4 Gt for the MFT, and 3.2 and 3.6 Gt for the COST-G and GSFC trends.Extending the integration area to a 50 km along the river yields a mass change for the MFT, COST-G and GSFC products of 53.5, 12.2 and 13.7 Gt, respecti vel y.This leads to a similar conclusion as in the Lake Victoria example.Extending the integration region to restore signal magnitude works because there are no other large geophysical mass change signals in the vicinity of the reservoir, but the superior spatial localization of the MFT is evident.The derived mass signal from the MFT closely matches the official reserv oir capacity, w hile the other two trend products underestimate the mass gain.Note that the ne gativ e ring around the reservoir is likely caused by ringing effects due to the limited resolution in the initial spherical harmonic parametrization.

High plains aquifer
The example shows North America's High Plains Aquifer (Ogallala Aquifer).The souther n par t of the aquifer showed a water mass loss in the last years, as can be seen both from the GRA CE/GRA CE-FO TWS observations and in situ groundwater well levels (see Fig. 7 ), provided by USGS (U.S. Geological Surv e y 2021 ) and analysed for the same time period as the GRA CE/GRA CE-FO trends.While TWS and groundw ater le vel data show dif ferent quantities and thus cannot be compared in absolute values unless converted to the same unit by applying information on the aquifer storage coefficient or specific yield, both data sets show similar spatial patterns of the trends, with a drying Southern part and a slightly wetting northern part of the aquifer.
The MFT data constrain the ne gativ e storage trend in the southern part of the aquifer much closer to its area of major groundwater depletion according to the in situ observations than the TWS trends based on the monthly GRA CE/GRA CE-FO records.

D I S C U S S I O N A N D C O N C L U S I O N S
Satellite g ravimetr y, as realized with GRA CE and GRA CE-FO , obtains information about the time-variable gravity field of the Earth from precise tracking data of two twin-satellites acting as a freefalling proof mass in a low orbit around the Earth.This implies a strong dependency of the observation accuracy on the spatial scales so that large scale signals are much better constrained than smaller spatial scales.In order to suppress spatially correlated errors at those smaller scales, either spatial constraints are used directly during the inversion of Mascon solutions, or spatial smoothing filters are applied to otherwise unconstrained spherical harmonics solutions.We have shown in this paper that both approaches are mathematically equi v alent (K usche 2007 ).
We also demonstrated in this work that the spatial resolution of gravity products derived from GRACE and GRACE-FO can be increased at the expense of temporal resolution.Through a new trend map based on an inversion of all available GRA CE/GRA CE-FO sensor data (mean field trend -MFT), we demonstrate that mass trend signals for several locations in the world with well-defined geographic boundaries (such as lakes, reservoirs, or regions with large groundw ater withdraw als for irrigation) are much better localized in this new MFT trend product when compared to trend pattern estimated a posteriori from monthly sampled GRACE products typicall y av ailable for geophysical applications.
In particular, we can separate the surface water storage trend observed with GRA CE/GRA CE-FO over Lake Victoria from the water storage trends of the surrounding area and lakes.The lake storage trend from satellite g ravimetr y of Lake Victoria matches  very well an independent estimate based on satellite altimetry.The benefit of the higher spatial resolution of the MFT product is also apparent for the water storage trends in Northwest India.In contrast to the standard monthly GRACE data sets, the new product can localize these trends to regions with pre v alent groundw ater irrigation.That supports the widely stated assumption that groundwater overuse in the region causes falling w ater le vels and, thus, a negati ve TWS trend.Onl y with the ne w MFT product can we accurately estimate the water storage gain of the ne wl y filled (2012-2016) Bo guchan y Reservoir in Siberia.Also for the High Plains aquifer, the higher spatial resolution of MFT allows for better localization of the trend signal inside the aquifer compared to conventional TWS products.The MFT approach can also be beneficial for groundwater data products as provided by the Global Gravity-based Groundwater Product (G üntner et al. 2023 ), where satellite g ravimetr y is combined with other remote sensing techniques and in situ data to separate water storage compartments.The spatial resolution of the GRA CE/GRA CE-FO contribution in this product is by far the lowest and could be brought closer to the other remote sensing data sets if only groundwater trends are of interest.
While MFT maximizes the spatial resolution at the expense of temporal information, we note that it is also possible to increase the temporal resolution at the cost of spatial resolution.An example of this are daily sampled estimates of the gravity field (Kvas et al. 2019 ).Here, additional stochastic constraints from independent geophysical models are required to stabilize the solutions as a consequence of the sparse satellite ground track coverage within a da y.How e ver, the dail y series published alongside the nominal monthly GRA CE/GRA CE-FO gravity fields from ITSG-Grace2018 have been successfully applied to identify major flood events (Gouweleeuw et al. 2018 ) and for the evaluation of atmospheric net-water fluxes (Eicker et al. 2020 ).Depending on the spatio-temporal characteristics of a mass transport process to be studied with satellite g ravimetr y, it is possible to optimize the results by trading spatial for temporal resolution and vice versa.A good example could be the estimation of coseismic signals, where the exact timing of the rupture is readily available from seismometers.In contrast, the spatial distribution of the mass shifts associated with the earthquake could be assessed from the difference in inverting all data before and after the event.Such an inversion was performed in the GOCO06s gravity field model (Kvas et al. 2021 ), where co-seismic gravity changes for three megathrust Earthquakes (the 2004 Indian Ocean earthquake, the 2010 Chile earthquake and the 2011 Tohoku-Oki earthquake) were co-estimated.Thus, tailoring the trade-off between spatial and temporal resolution to the signal of interest can bring new insights from satellite g ravimetr y beyond what standard monthly solutions offer.

A C K N O W L E D G M E N T S
This study received funding from the European Union's Horizon 2020 research and innov ation pro gramme under grant agreement no.870353, project G3P (Global Gravity-based Groundwater Product).This article was published Open Access with financial support from the University of Graz.
The CRediT Contributor Roles for this paper are as follows: AK: conceptualization, methodology, formal analysis, software, writing -original draft; EB: formal anal ysis, softw are, writing -original draft; HD: methodology, writing -original draft; AE: methodology, writing -original draft; TM-G: software; AG: methodology, funding acquisition, writing -original draft.

D ATA AVA I L A B I L I T Y
All data sets used in this study are publicly available, with data citations and repository locations referenced in the paper body.The MFT TWS grid is available under https://ftp.tugraz.at/outgoing/ITSG/MFT/ .

Figure 1 .
Figure 1.Comparison of linear TWS trend estimates based on monthly L3 products for the time span 2002-04 to 2021-01.

Figure 2 .
Figure 2. Difference between trend estimates based on Mascon and spherical harmonic parametrizations in the time span 2002-04 to 2021-01.

Figure 3 .
Figure 3. TWS trend estimate based on accumulated normal equations (mean field trend, MFT) for the period 2002-04 to 2021-01 (geocentre motion and postglacial rebound were corrected).Ocean bathymetry is adapted from NASA Visible Earth.

Figure 6 .
Figure 6.Krasnojarsk region in Siberia around Bo guchan y Reservoir: TWS trend from MFT, COST-G and GSFC Mascons for the period 2002-04-2021-01, dam location, reservoir inlet and reservoir extent.The Angara river along which the reservoir is located is highlighted in blue.

Figure 7 .
Figure 7. High Plains Aquifer in North America: TWS trend from MFT, COST-G and GSFC Mascons for the period 04/2002-12/2021.Circles give trends of in situ groundwater levels for several observation wells.Please note the different colourbars for TWS trends and groundwater level trends.