A hierarchical Constrained Bayesian (ConBay) approach to jointly estimate water storage and Post-Glacial Rebound from GRACE(-FO) and GNSS data

ABSTRACT Gravity Recovery and Climate Experiment (GRACE) and its Follow-On mission (GRACE-FO) have become an indispensable tool in monitoring global mass variations. However, separating GRACE(-FO) signals into its individual Terrestrial Water Storage Changes (TWSC) and surface deformation contributors, i.e. Post-Glacial Rebound (PGR), is desirable for many hydro-climatic and geophysical applications. In this study, a hierarchical constrained Bayesian (ConBay) approach is formulated to apply GRACE(-FO) fields and the uplift rate measurements from the Global Navigation Satellite System (GNSS) stations to simultaneously estimate the contribution of TWSC and PGR. The proposed approach is formulated based on a hierarchical Markov Chain Monte Carlo optimisation algorithm within a dynamic multivariate state-space model, while accounting for the uncertainties of a priori information and observations. The numerical implementation is demonstrated over the Great Lakes area, covering 2003–2017, where the W3RA water balance and the ICE-5G(VM2) and ICE-6G-D(VM5a) GIA models are merged with GRACE and GNSS data. Validations are performed against independent measurements, which indicate that the average root-mean-squares-of-differences between the PGR estimates and independent measurements reduced by after merging observations with models through ConBay. The ConBay updates, introduced to the long-term trends, as well as the seasonal and inter-annual components, are found to be realistic.

Therefore, reliable quantification of the hydrological cycle, its key fluxes and stores and its spatiotemporal variability is crucial for many applications related to hydrometeorology, water resources, and landatmosphere interactions, including flood and drought monitoring and prediction (Forootan et al., 2019Houborg et al., 2012;Li et al., 2019;Long et al., 2013;Slater et al., 2015), assessing water resources sustainability (Castellazzi et al., 2016;Forootan et al., 2014b;Scanlon et al., 2012), and identifying ecohydrological links between climate and vegetation (Singer et al., 2014).
Various hydrological models and hydrometeorological (drought/flood) monitoring systems exist to understand the temporal and spatial changes of water storage and water fluxes (Sood & Smakhtin, 2015). Even though the models aim for an adequate representation of the real world, uncertainties exist due to insufficient model realism, imperfect climate input data, inadequate empirical model parameters, and unrepresented feedbacks between model components (Rodell et al., 2004).
Today, organisations like the European Space Agency (ESA 1 ) and the National Aeronautics and Space Administration (NASA 2 ) provide various Earth Observation (EO) data that complement in-situ monitoring networks and represent a great potential to track changes in global water, carbon, and energy cycles (Lindersson et al., 2020). Among available satellite geodetic and remote sensing EO techniques, the Gravity Recovery And Climate Experiment (GRACE, 2002(GRACE, -2017 satellite mission (Tapley et al., 2004a,b) and its Follow-On mission (GRACE-FO, 2018 -onward) provide timevariable Earth's gravity fields that contain signals related to different processes such as non-steric sea level changes, Terrestrial Water Storage Changes (TWSC, i.e. a vertical summation of changes in water storage within plant canopies, surface water, snow, soil, and groundwater), ice sheet melting, and Post-Glacial Rebound (PGR), with a spatial resolution of a few 100 km and temporal resolution of , 10 days to 1 month in satelliteonly solutions (Flechtner et al., 2016).
The ability of GRACE(-FO) satellites to detect mass changes in the surface and sub-surface, which cannot be measured by any other satellite mission, and its sensitivity to water storage changes throughout all seasons, provides a unique opportunity to extract possible intensification's of the water cycle Kusche et al., 2016) and for drought monitoring at global (Forootan et al., 2019;Zhao et al., 2017b) and regional scales (Houborg et al., 2012;Schumacher et al., 2018a;Sinha et al., 2017;Thomas et al., 2014;Zhao et al., 2017a). Although GRACE(-FO) time-variable fields represent an accurate superposition of water storage changes and large-scale deformation signals, separating this integrated signal into its contributors is desirable for many geodynamic and hydro-climatic applications.
For the hydrological applications of previous signal separation studies, the effect of PGR is often removed from GRACE(-FO) data, during the post-processing, as a linear trend, using the output of a GIA model (e.g. Sasgen et al., 2012;Wahr & Zhong, 2012). However, large uncertainties in the present-day global GIA models exist due to the insufficient input data, e.g. uncertainty of ice loading history (Spada et al., 2011). This is also shown by Guo et al. (2012), who found considerable regional differences between 14 forward model solutions, disagreeing even on the sign of vertical land motion in some areas. Therefore, using GIA model outputs to remove the effect of PGR from GRACE(-FO) data may negatively affect the accuracy of hydrological mass estimations in signal separation approaches, for example, in regions such as North America, Greenland, and Scandinavia.
The loading deformation, due to PGR, is the result of the Earth's response for the redistribution of surface mass (Chanard et al., 2014;Craig et al., 2017;Fu & Freymueller, 2012), which can be observed by GNSS stations (D. Dong et al., 2002). Many studies have demonstrated that GNSS and GRACE provide comparable signals and the load effect derived from GRACE based on the layered elastic earth model can be a firstorder approximation of the seasonal displacement observed by GNSS (Argus et al., 2021;Chanard et al., 2018;Davis et al., 2004;Kusche & Schrama, 2005;Rietbroek et al., 2012;Tregoning et al., 2009;Van Dam et al., 2007;Wang et al., 2013). In the light of these investigation, few studies made use of both GRACE and GNSS to study seasonal water cycles (Argus et al., 2014a(Argus et al., , 2020Fu et al., 2015;Fu & Freymueller, 2012;Hao et al., 2016;Larochelle et al., 2018;Pan et al., 2016). These studies, however, did not address the separation of GRACE signal to its individual surface and sub-surface water storage components.
In this study, we propose a new hierarchical 'Constrained Bayesian (ConBay)' optimisation approach for a joint estimation of the land hydrology compartments and PGR rates from GRACE(-FO) field estimates, while in-situ GNSS measurements are used to constrain the estimated PGR rates from GRACE(-FO) data. In the ConBay approach, instead of removing PGR from GRACE (-FO) field estimates, during the post-processing, the outputs of a GIA model are used as a priori information, along with the hydrological model outputs, to simultaneously separate land hydrology and PGR from GRACE(-FO) data. In many glacial regions of the world, PGR signals are 'contaminated' by vertical elastic crustal deformation, which is induced by present ice mass change. Therefore, improvement in the estimation of PGR enhances the estimation of the magnitude and pattern of the elastic crustal deformation (surface deformation) within the continental regions.
It is worth mentioning here that the Ensemble Kalman Filter (EnKF) used for DA and C/DA (Schumacher, 2016) can also be classified as a Bayesian-based technique because the cost function for updating the conditionality of unknown state parameters on the measurement data is formulated based on the Bayes theory (see, e.g. Evensen, 2003;Fang et al., 2018;Schumacher, 2016). Schumacher (2016) discussed that the statistical information used in EnKF-DA is restricted to the covariance matrices of the observations and models. Therefore, the observation error model and the spatial resolution of GRACE TWSC have a significant influence on C/DA results. Instead of limiting the statistical information in the data to the use of their covariances, a sampling of their Probability Density Functions (PDFs) and a Bayesian optimisation approach is adopted in ConBay that results in more realistic estimations of states and their errors. More important, previous DA techniques rely on model equations to relate water and energy fluxes with water storage changes, and therefore, combining information from observations (e.g. GRACE(-FO) TWSC) and models is performed in a physically justified way. Unlike the previous DA, the proposed Bayesian technique in this study is implemented in an offline mode, where we do not need to run the models, and we only use the outputs of the available hydrological models to merge with the observations. The proposed ConBay approach is flexible to account for the uncertainties of observation and models and can be applied both globally and regionally and for different spatial and temporal resolutions, which does not have any influence on final results.
Particle Filter (PF) and Particle Smoother (PS) are other types of Bayesian approaches (Särkkä, 2013), which have been used in some hydrological applications such as Plaza Guingla et al. (2013); Weerts and El Serafy (2006) to assimilate the observations into the models. These techniques use a set of particles (also called samples) to represent the posterior distribution of some stochastic process given noisy and/or partial observations. However, since the computational cost of PF grows with the number of particles, choosing a specific number of particles in the design of filters is a key parameter for these methods. Therefore, PF and PS might not be efficient for high-dimensional fusion tasks, such as global hydrological application (e.g. Bain & Crisan, 2008;Snyder et al., 2008), while the proposed Bayesian approach in this study provides the ability to deal with high-dimensional fusion tasks, and its computational load is much lower than PF and PS.
In this study, to evaluate the performance of the proposed approach its implementation is presented within the Great Lakes area, the Unites States (US), during 2003 to 2017. The Great Lakes region that is located in the northeast of the US is selected as our case study, where mass changes due to PGR and surface deformations play dominant contributions in GRACE(-FO) data. Schumacher et al. (2018b) showed that large uncertainties exist between GIA models over this region making it a good candidate to study the performance of ConBay. Moreover, the groundwater storage is a major natural resource and links the Great Lakes and watersheds (Hall, 2006), whose interactions are not well understood (Sophocleous, 2002;Winter et al., 1998). This linkage needs to be better quantified before society could address some of the important water-resources issues in the region. Therefore, separating GRACE(-FO) signals into its hydrological and surface deformation compartments is of particular importance.
For this study, monthly GRACE field estimates, on 0:5 � � 0:5 � spatial grid points, and the in-situ GNSS vertical load rates (Schumacher et al., 2018b) are used as observations, and a priori information of hydrological changes is taken from the W3RA water balance model (Van Dijk, 2010) and PGR rates are taken from the ICE-5G (VM2) GIA model (Peltier, 2004). To evaluate the flexibility of the proposed approach to select various a priori information and to see how ConBay is sensitive to the use of different GIA models, we re-implemented the ConBay to merge GRACE and GNSS using ICE-6G-D(VM5a; Argus et al., 2014b;Argus et al., 2021;W. Peltier et al., 2018;W. R. Peltier et al., 2015) as a priori information of PGR rates. The ConBay-derived PGR rates and groundwater storage of both implementations (using ICE-5G(VM2) and ICE-6G-D(VM5a)) are then compared against independent in-situ USGS groundwater level observations and independent GNSS measurements (validation data sets), which were not used in the optimisation procedure.
To assess the impact of using GNSS to constrain PGR on the signal separation results, ConBay is reimplemented to merge only GRACE with the hydrological and GIA model outputs, while GNSS data were not used in the hierarchical level to accept/reject the estimated values of PGR. Therefore, the new implementation of ConBay, which is called here 'Unconstrained Bayesian-DA', is formulated similar to the MCMC-DA (Mehrnegar et al., 2020b). The PGR rates and groundwater storage changes derived from MCMC-DA (Unconstrained Bayesian-DA) and ConBay are then validated against independent validation data sets.

GRACE fields
In this study, we used the release six (RL06) monthly GRACE level 2 (L2) products, provided by the Center for Space Research (CSR 3 ) covering January 2003 -December 2017. The period is restricted to the availability of GNSS data from Schumacher et al. (2018b). The potential coefficients are truncated at the spherical harmonics of degree and order 90. In order to generate monthly Equivalent Water Heights (EWHs that contains both TWSC and PGR) from these products, recommended corrections are applied. Whereas, the degree 1 coefficients, which are not observed by GRACE, are replaced by those from Swenson et al. (2008) to account for the movement of the mass centre of the Earth. Degree 2 and order 0 (C 20 ) coefficients are replaced by more reliable estimates of the Satellite Laser Ranging (SLR) solutions following Chen et al. (2007).
The formulation in Wahr et al. (1998) is used to convert the L2 potential coefficients to 0:5 � � 0:5 � girded EWH fields within the Great Lakes (75 � W À 92:5 � W, 40 � N À 50 � N), while covering the period of 2003-2017. Correlated errors of the potential coefficients, which are caused by anisotropic spatial sampling of the mission, instrument noise, and temporal aliasing from incomplete reduction of shortterm mass variations (Forootan et al., 2014a), are reduced by applying the DDK2 filter (Kusche et al., 2009), which is comparable to a Gaussian filter of 340 km. The DDK filter is preferred here to the other filtering techniques since the final smoothed solutions are generally in better agreement with the TWSC output of global hydrological models (see, e.g. Werth et al., 2009). It also considers correlations between potential coefficients in a more rigorous manner compared to other filter techniques such as that of Swenson and Wahr (2006) who only model the order-dependent correlations.
For this study, instead of removing PGR, they are assumed to be unknown for the ConBay approach. Uncertainties of EWHs are computed by implementing a collocation error estimation Ferreira et al., 2016) using EWHs estimates from the CSR, Jet Propulsion Laboratory (JPL), and GeoForschungsZentrum (GFZ) L2 data.

W3RA water balance model
The Worldwide Water Resources Assessment (W3RA, Van Dijk, 2010)'s monthly averaged model states (snow, surface water storage, surface soil water (top layer), shallow-rooted soil water, deep-rooted soil water storage, and groundwater storage) are used as a priori information of TWSC components. For this study, the original code 4 is modified for the Great Lakes area by using daily 0:125 � � 0:125 � interpolated ERA-Interim reanalysis fields (Dee et al., 2011) of precipitation, albedo, 2-metre wind, as well as minimum and maximum temperature 5 as forcing data to run the model from 1980-2017. In W3RA, each cell is modelled independently of its neighbours, but lateral mass exchanges are accounted for by implementing a routing scheme. More details on the W3RA model can be found in Van Dijk (2010). Considering the resolution of GRACE, the original W3RA model outputs (0:125 � � 0:125 � ) are averaged on 0:5 � � 0:5 � grids, for the period January 2003 through December 2017.
Model uncertainty is estimated following Renzullo et al. (2014) by using the perturbed meteorological forcing approach. To this end, an additive error is assumed for the short-wave radiation perturbation of 50 Wm 2 , a Gaussian multiplicative error of 30% for rainfall perturbation, and a Gaussian additive error of 2 °C as the magnitude of the additive error air temperature perturbations. Estimated model uncertainty is used subsequently as the initial value of the variance/ covariance matrix of the unknown state parameter (see, Section 3.2.1) in the Bayesian inference, which is then updated through a forward-filtering and backward smoothing algorithm presented in Section 3.2.2.
Our motivation to select W3RA is its simplicity, which makes its computational load manageable for scientific applications (see examples of W3RA's applications in, e.g. Khaki et al. (2017); Forootan et al. (2019)), and its acceptable performance when compared with other commonly used Global haydrological or land surface models Schellekens et al. (2017).

Global Navigation Satellite System (GNSS) station data
In this study, in-situ GNSS vertical load rates released by Schumacher et al. (2018b) are used to estimate PGR rates in combination with GRACE data. This data set consists of 4072 in-situ stations over the globe (selected based on prior information from the GIA forward models to exclude tectonic signals), where 343 sites are located within the Great Lakes area, as our case study. The data covers the period of 2003-2017. Schumacher et al. (2018b) used the post-processed time series of the Nevada Geodetic Laboratory (NGL) as the starting point for providing an observational estimate of the GIA-related Vertical Land Motion (VLM). A fully-automatic post-processing strategy was then implemented to deal with outliers and jumps, atmospheric mass loading correction, elastic signal correction and filtering for stations where other sources of VLM are likely to dominate over GIA. To accurately account for the elastic response of the Earth's crust over Antarctica and Greenland, they used already corrected estimates that have dealt with the contemporary ice mass loading impact on elastic deformation, using high-resolution ice mass balance time-series (Khan et al., 2016;Martin-Espanol & Bamber, 2016). Moreover, they performed a spatial filtering strategy to select stations that are predominantly influenced by the long wavelength GIA signal and to exclude stations that are affected by local to regional hydrology (such as groundwater pumping).
To make the best of these in-situ measurements and to match the spatial resolution of GNSS data with GRACE EWHs, we applied the Least Squares Collocation (LSC) interpolation approach (Moritz, 1978) using the point-wise GNSS uplift rates (Figure 1(a)) and their uncertainties (Figure 1(b)) to obtain PGR uplift rates on a grid with 0.5-degree spatial resolution within North America (Figure 1(c)). The interpolated values over the Great Lakes area are used in this study, see, Figure 1(d). It is worth noting here that before implementing LSC technique to interpolate GNSS data on a grid, we excluded the observed values for 70 sites from 343 sites within the Great Lakes area, see, Figure 1(e). These stations are used as an independent validation to evaluate the ConBay-derived PGR estimates.
To compare gridded GNSS data with GRACE EWHs, PGR rates (uðθ; λÞ) are expanded to a series of spherical harmonics coefficients as where c u nm and s u nm are the associated cosine and sine coefficients, respectively, which can be obtained by solving a numerical integration (Wang et al., 2006) or a least squares approach (Sneeuw, 1994). Here, R is the mean equatorial radius, θ and λ are respectively the colatitude and longitude, and P nm represents the fullynormalised associated Legendre polynomials. The formulation in Wahr et al. (1998) is then applied to use these coefficients and to compute the EWHs that correspond to the PGR rates.

Global Isostatic Adjustment (GIA) model
The ICE-5G(VM2) global ice sheet reconstruction model output (Peltier, 2004), and -ICE-6G-D(VM5a; Argus et al., 2014b;W. Peltier et al., 2018;W. R. Peltier et al., 2015), which is one of the most recently published models of the GIA process in the ICE-NG(VMX) sequence from the University of Toronto, are used in this study as the a priori information for the PGR uplift rates within the Great Lakes area.
The ICE-5G(VM2) is a refined model of the Post-Last Glacial Maximum (Post-LGM) Global de-glaciation process, that has resulted through correction of the flaws in the ICE-4G (Peltier, 1994(Peltier, , 1996. ICE-5G differs significantly from the previous version (ICE-4G) at all Northern Hemisphere locations that were glaciated at Last glacial Maximum. These locations include all of northwestern Europe/Eurasia, the British Isles, Greenland and the North American continent (Peltier, 2004).
The ICE-6G-D(VM5a) model of the GIA process differs from all previous models in the ICE-NG(VMX) sequence in that it has been refined to constrain the thickness of local ice cover and the timing of its removal using the available Global Positioning System (GPS) measurements of vertical crustal motion in both Northern and Southern Hemispheres. Additional space geodetic constraints have also been applied for this model to specify the reference frame within which the GPS data are described (W. Peltier et al., 2018;W. R. Peltier et al., 2015).
The EWH estimates that correspond to the modelderived PGR rates are estimated following the procedure explained in Eq. (1; see, Figure 2(a,b)). The Root Mean Square of Differences (RMSD) between the mean of 14 GIA models in Guo et al. (2012) and the outputs of ICE-5G(VM2) and ICE-6G-D(VM5a) are considered as their uncertainties for the application part of this study. Considerable differences (up to 20 mm/yr in terms of EWHs rate) can be seen between ICE-5G(VM2) and ICE-6G-D(VM5a) within the Great Lakes area (Figure 2 (c-a,2-b)), specifically in North of the region. As it is expected, the ICE-6G-D(VM5a) GIA model shows more agreement (72%) with GNSS measurements, compared to the ICE-5G(VM2) GIA model. The reason is that ICE-6G-D(VM5a) reconciles virtually all available GPS measurements of vertical motion of the crust (W. Peltier et al., 2018;W. R. Peltier et al., 2015). Differences between these two models motivated us to use both GIA models for the application part of ConBay approach, to assess the performance of the proposed approach using various a priori information.

In-situ USGS groundwater level data
In this study, we use the data provided by the US Geological Survey (USGS) groundwater network, 6 which contains a record of groundwater levels between 1970-now for , 100,000 wells across the Conterminous United States (CONUS), which can be used as an independent validation data set to evaluate groundwater storage changes derived from W3RA model and/or ConBay results. The point-wise data for 2003-2017 are downloaded and filtered to exclude measurements with large data gaps (temporal gaps > 2 years), and those without meaningful variations, i.e. those time series which only contain linear and/or non-linear trends, without any other oscillations are excluded from data sets. Selected groundwater levels (, 763 wells mostly in southeast and southwest of the Great Lakes area) are then temporally averaged to produce monthly time series (see, Figure 3).
In-situ groundwater level data can be converted to groundwater storage changes using an effective Storage coefficient (Sc), where groundwater storage ¼ Sc � groundwater level. The USGS groundwater network covers a range of unconfined to confined conditions that are important to consider in evaluating groundwater level records and comparing with modelled groundwater storage changes. The storage coefficient (Sc), required to convert groundwater level to groundwater storage, can vary over several orders of magnitude from unconfined aquifers (e.g. between 0.02 and 0.3) to confined aquifers (known as storativity ,0:001; Freeze & Cherry, 1979). In systems with vertically stacked aquifers, it is often difficult to determine whether wells are screened in unconfined or confined aquifers, or both, which increases the uncertainty of estimating groundwater storage from groundwater level data. Therefore, different approaches have been introduced to approximate it regionally. For example, Rodell et al. (2007) applied an average value of 0.15 in the Mississippi River basin, and Xiao et al. (2015) used a range of 0.02 to 0.6 in the Mid-Atlantic Region of the CONUS based on the technical insights provided by USGS.
In this study, an approximate value of Sc (between 0.05 and 0.34) is estimated for each well based on the STATSGO soil texture class for the location of observed data. 7 The presented values of Sc for different soil texture class are extracted from Richey et al. (2015). These values are then used to convert USGS groundwater level to groundwater storage changes, which can be used as an independent validation data set to validate groundwater storage changes of ConBay. An overview of observation and models used in the application part of this study is presented in Table 1.

Method
ConBay, as an extension of MCMC-DA (Mehrnegar et al., 2020b), is formulated based on hierarchical multivariate state-space models (dynamic system), (I) between GRACE(-FO) observations and model outputs (both hydrological and GIA model), and (II) between GNSS measurements and the updated values of PGR derived from GRACE data in (I).
Learning dynamical systems (Thelen & Smith, 1998), also known as system identification or time series modelling, aim to create a model or improve an existent model based on measured signals (Lennart, 1999). In this study, the state-space models, also known as hidden Markov (Rabiner, 1989) and latent process models (Koller & Friedman, 2009), are used that describe the probabilistic dependence between the latent (unobserved) state variables (e.g. hydrological/GIA model outputs) and the observed measurement (e.g. GRACE EWHs). The Markov model is a pseudo-randomly stochastic process, which assumes that future states depend only on the current state, not on the events that occurred before it (it is known as the Markov property). This assumption enables reasoning and computation with the model that would otherwise be intractable (Gagniuc, 2017). The presence of a latent state, within the state-space model, allows for a succinct  representation of the dynamics in the form of a Markov chain. The state contains the information about the dynamic system (e.g. hydrological process within the Earth system), which is essential to determine future forecasts. Estimation of the unknown state parameters and the temporal dependency between them to solve the multivariate state-space model between GRACE (-FO) and model outputs are realised through a combination of a forward-filtering backwardsmoothing recursion approach (Kitagawa, 1987), and a Gibbs sampling (Gelfand & Smith, 1990;Smith & Roberts, 1993) algorithm. The central hypothesis for these formulations is that the magnitude of changes in water storage components depends on the history of hydrological processes. However, there is no or little physical knowledge about how this dependency varies over time. The combination of Gibbs sampling and forwardfiltering backward-smoothing approach allows a simultaneous dynamic estimation of the individual water storage compartments and temporal dependencies between them. Mehrnegar et al. (2020b) demonstrated that, compared to the DMDA approach (Mehrnegar et al., 2020a), which applies an empirically estimated constant value to control the temporal dependency in state-space model, a dynamic estimation of the unknown temporal dependency along with the unknown state parameters yields more realistic individual water storage estimates.
In a multivariate state-space model, when analysing multivariate data, one concerning issue is that the relationship between model parameters is unknown (or it is extremely difficult to determine). For instance, it could be expected that certain parameters have a larger effect on the dependent variables than other parameters. This issue in the application of separating TWSC and PGR from GRACE(-FO) estimates can be seen explicitly in the glacial regions, where PGR has a large effect on hydrological estimations. Moreover, PGR manifests as a trend in the relatively short era of the GRACE(-FO) mission, which needs to be introduced to the Bayesian separation framework. These theories can be transformed into the Bayesian fusion technique with specific inequality/equality constraints on the means and regression coefficients. Therefore, in the ConBay approach, a Metropolis-Hastings (Chib & Greenberg, 1995) algorithm is formulated, in a hierarchical level, to use in-situ GNSS measurements and constrain the GIA part by accepting or rejecting the updated value of the GIA model suggested by GRACE(-FO) data. The posterior estimated values of the state parameters are then used to update hydrological and GIA model outputs and their uncertainties.

Multivariate state-space model with unknown state equation
The objective of state-space modelling is to compute the optimal estimate of the state parameters given the observed data, which can be derived as a recursive form of the Bayesian rule (Brown et al., 1998;Z. Chen et al., 2010b). Multivariate state-space model between GRACE(-FO) observations and multiple a priori information (e.g. hydrological and GIA model outputs) can be represented by the observation and the state equations (Bernstein, 2005) as (2) In Eq.
(2), Y t ¼ ½y 1 ; y 2 ; . . . ; y P � t represents the vector of observation for P spatial grid points, at time t ¼ 1; 2; . . . ; T, while Z t and X t are two separate diagonal matrix to store a priori information from two different sources, and Θ t and β t are the unknown state parameters to make a relationship between the observation and a priori information.
For the application of ConBay in this study, Y t denotes the GRACE(-FO) field estimate, the diagonal elements of Z tðP�PÞ contain the hydrological model outputs, and the diagonal elements of X tðP�PÞ contains EWHs derived from the GIA model for the spatial grid point p ¼ 1; 2; . . . ; P at time t. Each of the diagonal elements of Z t is a 1 � K vector of ½z 1;p ; z 2;p ; . . . ; z K;p � t , where K is the number of individual water storage components, such as snow, canopy, surface water, soil water, and groundwater storage changes. Θ t is a P � 1 vector, where each element itself is a K � 1 vector containing the unknown state parameters for the water storage components, ½θ 1;p ; θ 2;p ; . . . ; θ K;p � T t , and β t is a P � 1 vector representing the unknown state parameters for the GIA signal, corresponding to the spatial grid points p ¼ 1; 2; . . . ; P.
A hierarchical constraint equation is formulated here to use the second observation data sets, i.e. rates from in-situ GNSS measurements, for controlling the sampling of β t derived from Eq. (2) as where G t is a P � 1 vector of the observation at spatial grid point p ¼ 1; 2; . . . ; P, and γ t represents their uncertainties. In Eq.
(2), (3), and (4), ε t , δ t , and γ t are additive innovations (i.e. residuals) corresponding to the observation equation, state equation, and the constraint equation, respectively. The distribution of ε t and γ t is assumed to be Gaussian with the mean value of zero and V t and U t to be the error covariance matrices, which vary over time. The state residual δ t is assumed stationary Gaussian distributed and independent from ε t and γ t , with the mean value of zero and an error covariance matrix of Q. Thus, the distribution of the additive innovations can be written as The uncertainty of the first set of observations (e.g. GRACE(-FO) measurements) is reflected in V t , while the uncertainty of the second set of observation (e.g. in-situ GNSS measurements) is reflected in U t . The error covariance matrix Q (corresponding to the state innovation δ t ) defines the temporal dependency between various compartments of our a priori information, which is unknown, and will be simultaneously estimated with the unknown state parameters Θ t and β t through the Gibbs sampling algorithm described in Section 3.2. The ConBay algorithm is formulated in the next section as a combination of a forward-filtering and backward-smoothing approach, Gibbs sampling, and Metropolis-Hastings algorithm to estimate the unknown state parameters Θ t and β t , and the covariance matrix Q, while the generated samples of β t in each iteration of Gibbs sampling are not accepted automatically as posterior samples; instead they are introduced as candidate samples to the hierarchical Metropolis-Hastings algorithm to be accepted or rejected based on the GNSS measurements. These candidate samples are accepted probabilistically based on the acceptance probability α (computed in Eq. (10)), which is estimated using the posterior probability distribution of β t conditional on the second set of the observations (e.g. GNSS measurement) based on the Eq. (4).

ConBay Algorithm
The multivariate state-space model, defined by Eq. (2) and (3), provides the conditional distribution of the parameter of interest, i.e. Θ t , β t and Q. To solve the state-space model, sampling techniques can be applied to generate samples from the posterior distribution of (i) time varying state parameters ðΘ 1:T ¼ ½Θ 1 ; Θ 2 ; . . . ; Θ T �Þ and ðβ 1:T ¼ ½β 1 ; β 2 ; . . . ; β T �Þ and (ii) the error covariance matrix of δ 1:T , i.e. Q, conditional on the first set of observation ðY 1:T Þ and its error covariance matrix (V 1:T ), and the rest of unknown parameters, i.e. Q in (i) as well as Θ 1:T and β 1:T in (ii).
Gibbs sampling (Gelfand & Smith, 1990;Smith & Roberts, 1993) is one of the most frequently used MCMC techniques to obtain samples from the posterior distribution. The idea in Gibbs sampling is to generate posterior samples by sweeping through each variable (or block of variables) to sample from its conditional distribution with the remaining variables fixed to their current values. For instance, consider the random variables ðX 1 ; X 2 ; . . . ; X D Þ.
Gibbs sampling is started by setting these variables to their initial values X This process continues until convergence (the sample values have the same distribution as if they were sampled from the proper posterior joint distribution).
Essentially, Gibbs sampling reduces the problem of sampling X :¼ ðX 1 ; X 2 ; . . . ; X D Þ to the problem of conditionally sampling of variables X d , for d ¼ 1; 2; . . . ; D. Since X d are of lower dimension (perhaps even one-dimensional), compared to X, they may be easier to sample by conventional methods. Therefore, when the conditional distribution of each variable is known, which is the case here, Gibbs sampling is an efficient approach to obtain samples from the posterior distribution.
However, generated samples from the posterior distribution of β t (e.g. state-space parameters corresponding to the GIA effects in this study) in step (i) are not accepted automatically as posterior samples and will be controlled by the second set of observations G t (e.g. insitu GNSS measurement) based on Eq. (4). Before generating samples of covariance matrix Q, a Metropolis-Hastings algorithm is formulated in a hierarchical level to estimate the acceptance probabilities using the posterior distribution of candidate samples β t conditional on the observations G t . The acceptance probabilities are then used to determine the best generated samples of β t , which leads to decrease the RMSD between the updated value of X t (i.e. X t β t ) and the observations G t according to Eq. (4). The accepted values of the generated samples β t are then used to generate the posterior samples of Q in step (ii).
The state-space model given by Eq. (2) and (3) is linear, and it is assumed that the distribution of the observations and the a priori information are Gaussian and independent from each other (though their covariances are assumed to be spatially correlated). Therefore, the conditional posterior distribution of Θ 1:T and β 1:T are the products of Gaussian Probability Density Functions (PDFs), which can be generated using a standard simulation smoother introduced by Carter and Kohn (1994). Samples generated from the conditional posterior of Q are the product of independent Inverse-Wishart distributions (Schuurman et al., 2016), which are defined on symmetric and positive definite matrices and used generally as the conjugate prior for the covariance matrix of a multivariate normal distribution in the Bayesian inference. The mathematical formulations to estimate the posterior distribution of the parameters of interest are explained in details in what follows.

Step1:
Define initial states, or prior values, for the unknown state parameters Θ t and β t , where t ¼ 0, and the prior value of the error covariance matrix of additive innovation δ t , i.e. Q ðiÞ , where i denotes the iteration number in the Gibbs sampling (it is zero here). Details of initial values for the unknown parameters are explained in Section 3.2.1, where the prior values for β t are determined similar to those of Θ t .
Step2:  In Eq. (7) and for the rest of the equations, pð�j�Þ is used to denote a generic PDF of (a variable such as � ) conditional on (another variable such as � ), while N (.) indicates a Gaussian PDF, and � (.) is an operator to multiply PDF. A Forward-filtering backward-smoothing approach, as in Kitagawa (1987), is used to estimate the unknown state parameters Φ tjtþ1 =[Θ tjtþ1 ; β tjtþ1 � and their error covariance matrices � tjtþ1 . Details and the corresponding equations are provided in Section 3.2.2. The outputs of the forward-filtering backwardsmoother approach are the generated samples of conditional on the GNSS measurements G t and its error covariance matrix U t as where X ðjÞ t is the updated value of X t (GIA model output in this study) and is estimated using the mean value of generated sample β ðjÞ t , while � ðjÞ X t denotes the uncertainty of X ðjÞ t , which is estimated using the error propagation formula as where j ¼ 1; 2; . . . ; i. It is worth noting here that The minimum value of α ðlÞ t indicates that the generated samples of β ðlÞ t , and therefore, the update value of priori information X t , i.e. X t β ðlÞ t , is more fitted to the observation G t (according to the Eq. (4)), compared to those of X t β ðiÞ t . β accepted ðiÞ 1:T derived from this step, along with the Θ ðiÞ 1:T of Step 2 are then used to generate samples of the unknown covariance matrix of additive innovations δ t , i.e. Q in Eq. (5), in the next step of the Gibbs sampling.
Step4: Sample Q ðiÞ from the posterior PDF of Q ðiÞ conditional on the observed data ðY 1:T Þ and its error covariance matrix (V 1:T ), and Φ where In Eq. (15), IW (.) denotes an Inverse-Wishart PDF, � Q is the posterior scale matrix, and ν is an initial value that is chosen as the degree of freedom to define the conjugate prior for Q as the product of independent Inverse-Wishart distribution (see, Section 3.2.1), and � ν is the posterior value of the degree of freedom. In all these equations, T denotes the total number of time steps t.
Step5: Return to step 2 and continue the iteration until a breaking criterion is satisfied. Two important issues that must be addressed while implementing MCMC algorithms are where to start and when to stop the algorithm. These two tasks are related to determining convergence of the underlying Markov chain to stationarity and convergence of Monte Carlo estimators to population quantities, respectively. It is known that under some standard conditions on the Markov chain, for any initial value, the distribution of unknown parameters converges to the stationary distribution (see, e.g. Meyn & Tweedie, 2012;Robert & Casella, 2013) as i ! 1.
Therefore, MCMC algorithms are typically run for a large number of iterations (in the hope that convergence to the target posterior will be achieved). In this study, a simple graphical method suggested by Brooks and Roberts (1998) and Sinharay, (2003) used to determine the number of iterations and to define the convergence of the sampling algorithm. This is done by creating a time-series plot for each of the parameters of interest, i.e. Θ t and β t in Eq. (2), and Q in Eq. (5), to view the path traversed by the chains. From the obtained results (figures not shown), it can be found that after 10,000 iterations, the Gibbs sampling converged to stationary distributions. However, to increase confidence in the process, it is required to select more iterations, i.e. N=20,000. Moreover, the early M=500 iterations are discarded as the 'burn-in' period. The reason for this is that after initialising the sampling algorithm with a priori values for the unknown state parameters Θ t and β t , and the covariance of additive innovations Q, samples from early iterations may not necessarily be representative of the actual posterior distributions.

Specifying Initial Values for Unknown Parameters
The first step of Gibbs sampling to solve the statespace model, given by Eq. (2) and (3), is to define initial states, or prior values for the unknown state parameters Θ t and β t , where t ¼ 0, and the error covariance matrix of additive innovation δ t , i.e. Q ðiÞ , where i denotes the number of iteration in the Gibbs sampling (it is zero here). As discussed in Section 3.2, since the state-space equation is linear, and with assuming that the distribution of observations and models, i.e. Y t , Z t , and X t , to be Gaussian and independent from each other, the conditional posterior distribution of Φ  (7) is a product of their (Gaussian) PDFs. Therefore, the initial value of Φ 0 is Gaussian distributed and is defined as where Nð:Þ represents a Gaussian (normal) distribution, and Φ 0j0 ¼ ½Θ 0j0 ; β 0j0 � and � 0j0 ¼ ½� Θ 0j0 ; 0; 0; � β 0j0 � are the means and variances of Θ 0 and β 0 . In GRACE(-FO) applications, Φ 0j0 is chosen to be 1, because in theory the summation of individual water storage values Z t and EWHs of PGR of X t must be equal to the GRACE(-FO) field estimates Y t in Eq. (2). Uncertainty of a priori information (e.g. the uncertainties of simulated individual water storage of W3RA (Section 2.2), and the GIA model output (Section 2.4) are stored in � 0j0 , which are needed for computing updates by the Gibbs sampling in Step 2. In Eq. (16), Q is the error covariance matrix of additive innovation δ t in Eq. (3) and defines the temporal dependency between water storage states Θ t at each time point to the previous time steps. Considering Eq. (15), the conjugate prior for Q can be estimated by an independent Inverse-Wishart distribution (IWð:Þ) as Q ð0Þ , IWðS; νÞ ; In Eq. (18), S is a K � K scale matrix, where K is the number of unknown state parameters, and ν is the degree of freedom Schuurman et al. (2016). S is used to position the Inverse-Wishart distribution in the parameter space, and ν > K þ 1 set the certainty about the prior information in the scale matrix. In this study, ν is set to be K þ 1, which is the minimum value that can be chosen for this parameter (Primiceri, 2005). Schuurman et al. (2016) compared three Inverse-Wishart prior specifications: (I) a prior specification that is based on an identity matrix, and is often used as an uninformative prior in practice, (II) a data-based prior that uses input from maximum likelihood estimations, and (III) the default conjugate prior proposed by (Kass & Natarajan, 2006). Their results showed that the data-based maximum likelihood prior specification for the covariance matrix of the random parameters, based on estimates of the variances from the data, performed the best, compared to the other techniques. They also found that when the prior is specified too far from zero (e.g. Inverse-Wishart prior with S as an identity matrix), this will result in an overestimation of the variances. However, specifying the central tendencies too close to zero will result in an underestimation of the variances, firstly because too much weight is shifted towards zero, and secondly because an element of the scale matrix set close to zero will also have a small variance. Following Cogley (2005), Primiceri (2005), and Schuurman et al. (2016), the scale matrix S is chosen to be a constant fraction of the variance of the initial values Θ 0 as K 2 Q � 0j0 . Therefore, Eq. (18) is rewritten as Q ð0Þ , IWðK 2 Q � 0j0 ; 1Þ : In Eq. (19), � 0j0 is the covariance matrix that is derived from the ensemble of W3RA (see, Eq. (17)). Following Primiceri (2005), K Q is chosen to be 0:01, which allows Θ t to be a temporal variable. More details about selecting different initial values and their impact on the Inverse-Wishart estimation can be found in Schuurman et al. (2016). Initial values chosen for Θ 0 and Q ð0Þ are then used in Eq. (7) and (16) within Step 2 and Step 4 of Gibbs sampling, respectively.
The forward-filtering backward-smoothing is a recursive approach that consists of two steps: (1) Forward-filtering: It consists of a standard Kalman filter procedure as in Carter and Kohn (1994) to recursively estimate Φ tjt and � tjt , for t = 1,2,.,T, given the initial values of Φ tjt and � tjt when t = 0 (Φ 0j0 and � 0j0 , see, Section 3.2.1). The Kalman filter recursion is presented in the following equations as The last elements of the Kalman filter recursion when t = T, i.e. Φ TjT and � TjT , are used to generate the samples of Φ T in Eq. (22) as where N(.) denotes the normal distribution of Θ T with the mean value of Θ TjT and the variance of � TjT .
(2) Backward-smoothing: The outputs of step (1) (i.e. Φ tjtÀ 1 , � tjtÀ 1 , Φ tjt , and � tjt for t = 1,2, . . .,T), and the generated sample of Φ T , derived from Eq. (22) are then used in Eq. (23) to update Φ tÀ 1jt and � tÀ 1jt , and generate samples of Φ tÀ 1 , for t ¼ T; T À 2; . . . ; 1. For a generic time t, the updating formulas of the backward recursion smoother can be written as The backward recursion smoother is started from time T and continues until time 1. The output of the forwardfiltering backward-smoother approach is the generated samples of Φ 1:T , i.e. Θ 1:T and β 1:T , conditional on the observed data ðY 1:T Þ and its error covariance matrix (V 1:T ), and the unknown covariance matrix of the additive innovation Q, which is defined as pðΦ ðiÞ t jΦ ðiÞ tþ1 ; Y 1:T ; V 1:T ; Q ðiÀ 1Þ Þ in Eq. (7). The generated samples of Θ 1:T and β 1:T are then used in the third and forth step of Gibbs sampling to generate samples of the unknown covariance matrix of additive innovations, i.e. shown by Q in Eq. (15).

Updating a priori information and their uncertainties by ConBay
At the end of the ConBay algorithm, N À M generated samples for Θ 1:T and β 1:T , are used to estimate the posterior value of unknown state parameters � Θ 1:T , and � β 1:T as where � Θ 1:T , and � β 1:T are then used to update the a priori information Z t (e.g. hydrological model outputs) and X t (e.g. GIA model output) as The uncertainties of the ConBay updated signals are estimated using the variance of generated samples, i.e.
� ðiÞ Θ t and � ðiÞ β t in Eq. (9), which is derived from the forward-filtering backward-smoothing approach following the equations in Section 3.2.2. To this aim, the posterior value of � � Θ t and � � β t for t ¼ 1; 2; . . . ; T are obtained as The diagonal elements of the error covariance matrix of � � Θ t contain δ � θ 2 k;t corresponding to each compartment of the a priori information Z t , k ¼ 1; ::; K, which can be used to estimate the uncertainties of Ẑ t through an error propagation procedure according to the Eq. (27). δẑ 2 k;p;t ¼ δ � θ 2 k;p;t :z 2 k;p;t þ δz 2 k;p;t : � θ 2 k;p;t ; where δẑ 2 k;p;t are the uncertainties of the ConBay updated values of a priori information (ẑ k;p;t ) for p ¼ 1; . . . ; P, t ¼ 1; 2; . . . ; T and k ¼ 1; 2; . . . ; K. The uncertainty of X t derived from Eq. (25) is estimated similar to that of Ẑ t based on Eq. (27). The work-flow of the ConBay approach is summarised in Figure 4.

Departure from the unconstrained Bayesian formulations
The MCMC-DA (Mehrnegar et al., 2020b) which is formulated based on the state-space model (Eq. (2), (3)) is a good candidate of an unconstrained Bayesian approach to separate GRACE(-FO) signals without considering extra information about PGR rates. In fact, MCMC-DA and ConBay are the same is the Steps 1, 2, 4, and 5 of Section 3.2 to estimate unknown state parameters and temporal dependency between them, while religiously account for the uncertainties of the observations and a priori information. Their main difference is that the hierarchical Metropolis-Hasting algorithm (see Step 3 in Section 3.2) to constrain the PGR rates based on GNSS measurements is not implemented in MCMC-DA. Thus, making MCMC-DA an unconstrained Bayesian technique. A comparison between these two techniques will show how beneficial is the joint application of GRACE and GNSS data in signal separation studies.

ConBay application within the great lakes area using ICE-5G(VM2) GIA Model
Linear rates of PGR in terms of EWHs are presented in Figure 5, where considerable differences (, � 14 mm/ yr) are found between the ICE-5G(VM2) GIA model outputs (Figure 5(a)) and GNSS measurements ( Figure  5(d)) in the northeast and the northwest. The model PGR rates are found to be under-estimated (compared with the observations in the northeast, and overestimated in the northwest of the Great Lakes (see, Figure 5 (d-a)). In the northeast, the rates of EWHs of PGR derived from GNSS measurements are estimated around 32 mm/yr, while in northwest these values are estimated up to 4 mm/yr. ICE-5G(VM2) GIA model, however, simulates these values to be between ,15 to 20 mm/yr.
The updated values of PGR derived from both MCMC-DA ( Figure 5(b)) and ConBay (Figure 5 (c)) show more agreements with GNSS data, compared to the original model outputs. Considerable improvements are seen in Figure 5 (d-b,5-c), where the average RMSD between GNSS measurements and GIA model output is reduced from 15 to 2 mm/yr in the northeast, and from À 17 to À 4 mm/yr in the northwest after implementing both MCMC-DA and ConBay.
Comparing MCMC-DA and ConBay in Figure 5 (d-b,5c) shows that using in-situ GNSS measurements to constrain the estimation of PGR can reduce the effect of uncertainties in the middle of Great Lakes, where differences between MCMC-DA and in-situ GNSS measurements are found to be 13 mm/yr, in terms of EWHs, while those between ConBay and in-situ measurements are estimated to be up to 5 mm. It is worth noting that in the middle of Great Lakes, shown by the grey patched area in Figure 5 (d-b,5-c), the estimated PGR rates, derived from both techniques, are over-estimated, compared to the in-situ GNSS measurements. This is likely due to the gaps in the GNSS data. Interpretation of the results within this area might be done with cautions.
To validate the PGR rates of MCMC-DA and ConBay, the GNSS measurements of 70 stations ( Figure 5(e)) are used, which were excluded from GNSS data before implementing the ConBay approach. Validation is done in terms of differential values between GNSS measurements and GIA model outputs before and after implementing MCMC-DA and ConBay (see,   5 (e-a,5-b,5-c), respectively). The averaged RMSD between the GIA model and GNSS measurements are estimated to be 7.8 mm/yr, while after implementing MCMC-DA and ConBay this value is reduced to be 4.8 and 2.3 mm/yr, respectively. The results indicate that both methods decrease the uncertainties but that of ConBay is closer to GNSS observations. In order to test whether the implemented mergers introduce unwanted anomalies to the a priori information, we estimate the Root Mean Square of Errors (RMSE) of estimated PGR components, derived from MCMC-DA and ConBay, after removing the long-term linear rates, see, Figures 6 (a,b), respectively. The results indicate that in fact unwanted signals with an amplitude of 30 À 50 mm can be seen in the PGR estimates of MCMC-DA. However, in ConBay, the constraint equation and GNSS data prevent this error and the magnitude of unwanted signal is found to be of the level of GRACE noise (less than ,15 mm). It shows that in an ideal situation, where the input observation and models are perfect, without any errors, the uncertainties of the ConBay results shall be close to zero. However, due to different reasons, such as imperfect model simulations (Rodell et al., 2004) and the errors of observations, it would be impossible to update a priori information without having any errors.
Before concentrating on the hydrology-related results of ConBay, in Figure 7, we provide an overview of the variance between GRACE EWHs (contain both hydrology and PGR) and modelled EWHs (as a summation of W3RA TWSC and EWHs of PGR derived from ICE-5G(VM2) GIA model) within the Great Lakes area between 2003-2017. Considerable differences are found between the measured and modelled EWHs, where the linear trends are shown in Figure 7 (a1,)) and the annual-amplitudes are presented in Figure 7 (a2,). For instance, the linear trends fitted to the GRACE EWHs in the northeast of the Great Lakes area are estimated to be ,32 mm/yr, which are much greater than those estimated by models (summation of W3RA TWSC and ICE-5G(VM2) EWHs of PGR) with the mean value of ,20 mm/yr. Further differences can be seen in the northwest, where GRACE shows positive trends up to ,10 mm/yr, while the summation of W3RA and GIA model shows ,18 mm/yr. In Figure 7 (c1), it can be seen that ConBay reduces the differences by 81% (from 11 mm/yr to 2 mm/yr on average) in the north part of the Great Lakes area.
Large differences in terms of annual-amplitude between measured (GRACE) and modelled EWHs can be seen in the northeast and southwest of Great Lakes, where those of GRACE are estimated up to 70 mm and those of W3RA are estimated to be greater than 90 mm in northeast and less than 50 mm in southwest. Differences in the seasonality of measured and modelled water storage changes, can be related to the errors in the forcing data and uncertainty in the model parameters to control these values (Van Dijk et al., 2011). The mean of annual-amplitudes fitted to the ConBay EWHs (Figure 7 (c2)) in northeast and southwest of Great Lakes are estimated to be ,65 mm, which is more close to those of GRACE (i.e. ,70 mm), compared to the original model outputs. To evaluate the performance of the ConBay to improve the estimation of EWHs, the RMSD between measured and modelled EWHs (Figure 7 (e1)) are compared with the RMSD between GRACE and ConBay EWHs (Figure 7 (e2)). The results indicate that the mean value of RMSD between measured and modelled EWHs is reduced by 63% (from ,68 to ,25 mm) within the Great Lakes area. Therefore, as expected, merging models with Figure 6. RMSE of EWHs of PGR, after remove PGR rates, derived from (a) MCMC-DA approach, i.e. without using GNSS measurements to constrain PGR estimates, and (b) ConBay approach that simultaneously uses GRACE and GNSS data for jointly estimating water storage and PGR rates. GRACE and GNSS measurements through ConBay approach improves the seasonal and inter-annual components of hydrological signals and reduces biases between the modelled and measured EWHs. Both MCMC-DA and ConBay are implemented in this study to simultaneously estimate land hydrology and surface deformation signals (PGR uplift rate) within the Great Lakes area. Therefore, it is expected that using GNSS measurements in ConBay to constrain the estimation of PGR rates affect the estimation of hydrological signals. To assess this hypothesis, the long-term linear trends and annual-amplitudes fitted to the groundwater storage changes derived from W3RA, MCMC-DA (un-constrained Bayesian approach) and ConBay are compared against in-situ USGS groundwater observation, as an independent validation data set in Figure 8. The RMSD and temporal correlation coefficients between in-situ USGS measurements and groundwater storage changes of W3RA, MCMC-DA and ConBay are used as two different metrics to validate our results. W3RA water balance model simulates groundwater depletion in southeast of Great Lakes with the rate of À 4mm=yr, on average, and groundwater recharge in northeast with the rate of 2:2mm=yr between 2003-2017 (Figure 8 (b1)).
After implementing both MCMC-DA and ConBay (see, Figures 8 (c1,)), the rate of changes in modelderived groundwater storage are reduced by ,52%, from À 3:8 to , À 1:8 mm/yr, on average, in southeast of Great Lakes. USGS groundwater storage also shows small negative trends in this region (between 0 and , À 1 mm/yr), which can be an evidence for the improved estimation of groundwater after merging model with GRACE and GNSS observations. The strong annual-amplitudes for groundwater storage changes are simulated by W3RA water balance model within the centre and south of the Great Lakes area (up to 60 mm, Figure 8 (b2)), which are reduced by ,33% after implementing MCMC-DA and ConBay (Figure 8 (c2,), respectively). The mean of RMSD between USGS observations and W3RA groundwater storage (Figure 8 (b3)) is estimated to be 44:21 mm. This value is reduced by 26% after merging W3RA and GIA model outputs with GRACE observations through MCMC-DA approach (Figure 8 (c3)). This improvement, in terms of the reduction in RMSD, is estimated to be 36%, when ConBay in implemented to merge GRACE and GNSS with the model outputs.
Positive correlation coefficients with the maximum value of ,0:9, mostly in southeast of Great Lakes, and the mean of 0:32 are estimated between USGS and W3RA groundwater storage in 95% of the validation points (723 wells from 763), which are not changed considerably after implementing the MCMC-DA (the mean value of 0.34) and ConBay (the mean value of 0.38) approaches.
Considerable differences between MCMC-DA and ConBay can be seen in the centre and southwest of Great Lakes, where a mixture of strong positive and negative trends ( � 6 mm/yr) can be seen in MCMC-DA groundwater storage. In these regions, we can also see a large RMSD, up to 60 mm, between MCMC-DA and in-situ USGS groundwater observations, which are considerably reduced (by 70%) when ConBay is implemented to separate land hydrology and PGR from GRACE and GNSS measurements. Therefore, it can be said that using in-situ GNSS measurements to constrain the PGR rates not only improves the estimation of PGR, but also influences the estimation of hydrological signals in terms of both linear trends and seasonal and interannual amplitudes.
The reason of finding differences between MCMC-DA and ConBay is the strong contribution of PGR in the GRACE signal and the known fact that these values are uncertain. Therefore, within the application of ConBay and MCMC-DA, to optimise updated value of the modelled EWHs at each time step, with respect to GRACE data, the unknown state parameters corresponding to the GIA model output (β t in Eq. (2)) gains the highest value of the posterior probability distribution compared to those of hydrological compartments (θ t in Eq. (2)). Therefore, the largest updates are introduced to the GIA model output in order to reduce RMSD between GRACE and model EWHs. As already discussed in Figure 6, without using in-situ GNSS measurement to constrain the estimation of PGR (as is the case in MCMC-DA), the PGR accepts a large portion of the update because statistically it will be more likely that the biases between a priori information and observations will be reduced. This update therefore does not care about the physical representation of PGR, which should be purely linear. As a result, the MCMC-DA estimation of PGR contains seasonal component and its hydrological estimation is negatively affected as the validation in Figure 8 represents. In order to visualise differences between the estimated parameters from MCMC-DA and ConBay, a comparison of posterior values of the unknown states θ soilwater t , θ groundwater t , and β PGR t are shown in Figure 9. The results correspond to only one spatial grid point with the latitude of 44 � and the longitude of À 90 � , for 2003-2017. These parameters (and others computed for other grids) can be used in Eq. (25) to update soil water storage, groundwater storage, and PGR rates, respectively. which can be seen in Figure 8. This shows that introducing GNSS observations in an extra step controls the updates and consequently improves the estimation.

ConBay application within the Great Lakes area using ICE-6G-D(VM5a) GIA model
To demonstrate how selecting different a priori information affects the final signal separation results, the ICE-5G(VM2) GIA model output is replaced by -D(VM5a; Argus et al., 2014b;W. Peltier et al., 2018;W. R. Peltier et al., 2015) for the re-implementation of ConBay over the Great Lakes area and its results are compared with those derived from first implementation, when we used the ICE-5G(VM2) GIA model (see, Section 4.1).   (a1) and (b1), and their differences with those of derived from ConBay using ICE-5G(VM2) are shown in (a2) and (b2), respectively. Differential values between GNSS measurements (validation points) and PGR rates derived from ConBay are shown in (a3). The values of RMSD between USGS and ConBay-derived GWS, in case of using ICE-6G-D(VM5a) GIA model, are shown in (b3).
The rates of new-derived EWHs of PGR and groundwater storage changes are shown in Figure 10 (a1,), respectively, where their differences with those derived from first implementation are shown in Figure 10 (,). The obtained results indicate that the differences between final signal separation results derived from both experiments are negligible (close to zero) for more than 90% of region.
Validations against in-situ GNSS measurements and USGS groundwater storage changes (see, Figure 10 (a3,), respectively), confirm that although ICE-6G-D (VM5a) shows more agreement with GNSS data compared to the ICE-5G(VM2; see, Figure 2), the RMSD between validation points and the ConBay-derived PGR and groundwater storage changes are not considerably changed for the new implementation (compare Figure 5 (e-c) and Figure 8 (d3)).

Summary and conclusion
In this study, a novel hierarchical Constrained Bayesian (ConBay) approach is formulated for a joint estimation of the land hydrology and surface deformation by merging the GRACE(-FO) and in-situ GNSS measurements, where hydrological and Global Isostatic Adjustment (GIA) models are used as a priori information.
ConBay is formulated as a combination of forwardfiltering backward-smoothing recursion approach (Kitagawa, 1987), the Gibbs sampling (Gelfand & Smith, 1990;Smith & Roberts, 1993), and the Metropolis-Hastings algorithm (Chib & Greenberg, 1995), where the first builds the relationship between observations and a priori information, the second accounts for unknown temporal dependencies between model states, and the latter implements the constraint equation. The proposed approach is formulated to account for the full error covariance matrix of the observed GRACE(-FO) data, the error estimates of GNSS data, and the uncertainty of the hydrological and GIA model outputs. Moreover, the proposed approach is flexible to implement both globally and regionally and for different spatial and temporal resolutions.
As a case study, the Great Lakes area and the period of 2003-2017 are chosen to merge GRACE and GNSS measurements, while the W3RA water balance model and the ICE-5G(VM2) GIA model outputs are used as the a priori information for hydrological and surface deformation signals. To assess the impact of the constraint equation on the signal separation results, MCMC-DA (unconstrained Bayesian-DA; Mehrnegar et al., 2020b) is performed to separate land hydrology and surface deformation from GRACE field estimates, where GNSS measurements are not used to constrain the estimation of PGR from GRACE data. Moreover, to see that how ConBay is sensitive to select different a priori information, the ConBay approach is reimplemented to merge GRACE and GNSS using ICE-6G-D(VM5a), instead of ICE-5G(VM2), as the a priori information for the PGR rates. The ConBay approach is tested by performing various comparisons between PGR rates, EWHs and groundwater storage derived from the original model outputs, ConBay, and MCMC-DA results. Validations are done against the in-situ US Geological Survey (USGS) groundwater storage observations and independent in-situ GNSS measurements which are not used in the optimisation procedure.
The numerical results indicate that after implementing ConBay the RMSD between model-derived EWHs and GRACE EWHs reduces by ,63%, on average, and the bias between GIA model output and the in-situ GNSS observation reduced by ,72% (from 7.8 to 2.2 mm/yr). MCMC-DA (unconstrained Bayesian-DA) results show the same improvement in the estimation of EWHs, but the bias between GIA model output and in-situ GNSS measurements reduced by 40% (from 7.8 to 4.7 mm/yr for the rate of equivalent water heights corresponding to PGR) after implementing MCMC-DA, which is considerably smaller than those of ConBay.
The estimated values of PGR through MCMC-DA contained an unwanted signal (see, Figure 6), which is not physically realistic. ConBay reduces this uncertainty by 67% (from 42 mm to 14 mm, on average), using the insitu GNSS measurements to constrain the updated value of the ICE-5G(VM2) GIA model.
Validation of groundwater estimates, derived from MCMC-DA and ConBay, against in-situ USGS groundwater observations provides an evidence for the hypothesis that after using in-situ GNSS measurements in the introduced Bayesian signal separation approach, the estimation of hydrological signals is improved in terms of both linear trend and seasonality. For example, we find that the RMSD between USGS observations and groundwater storage changes are reduced from 44.21 mm to 32.65 mm (26% improvement) after merging W3RA and GIA model with only GRACE data through MCMC-DA approach. This improvement, however, is greater (around 36%) when GNSS measurements are used in ConBay approach. From the results, it can be concluded that using GNSS constraint equation in ConBay defines an upper and lower boundary to update the GIA model output, which allows us to introduce more realist updates to the hydrological signals, compared to those of MCMC-DA. Therefore, using the GNSS constraint equation not only optimises the estimation of the land surface deformation, but also improves the estimation of the hydrological signals, in terms of both linear trends and seasonality.
The ConBay-derived PGR rates and groundwater storage changes of both implementations (using ICE-5G(VM2) and ICE-6G-D(VM5a)) are compared against validation data sets and the results indicate that the separated hydrological and PGR components from GRACE signal are not considerably changed in these two experiments. Therefore, we conclude that the ConBay approach is flexible to use different a priori information, and the Bayesian optimisation algorithm works well to introduce realistic updates to the background model. We also find that selecting ICE-6G-D (VM5a) leads to achieve faster convergences of Gibbs sampling (,5000 iterations), compared to those of ICE-5G(VM2) (,8000 iterations). Accordingly, it can be said that although the final signal separation results are independent of selecting the background models, the speed of Gibbs sampling convergence is influenced by this selection.