Variance–covariance analysis of two high-resolution regional least-squares quasi-geoid models

This paper investigates the full variance–covariance (VC) matrix of two high-resolution regional quasi-geoid models, utilizing a spherical radial basis function parameterization. Model parameters were estimated using weighted least-squares techniques and variance component estimation (VCE) for data weighting. The first model, known as the “RCR model,” is computed through the remove–compute–restore method, incorporating various local gravity and radar altimeter datasets. The second model, the “combined model,” includes the GOCO05s satellite-only global geopotential model as an additional dataset with a full-noise VC matrix. Validation of the noise VC matrix scaling for each quasi-geoid model is achieved by comparing observed and formal noise standard deviations of differences between geometric and gravimetric height anomalies at GPS height markers in the Netherlands. Analysis of the noise VC matrix of height anomalies at grid nodes reveals significantly smaller formal noise standard deviations for the RCR model compared to the combined model. This difference is attributed to VCE assigning larger weights to the GOCO05s dataset, which exhibits greater noise standard deviations for the specific spatial scales used. Additionally, the formal noise standard deviations of height anomaly differences, relevant for GNSS-heighting, favor the RCR model. However, the disparity between the two models is smaller than implied by the height anomaly noise standard deviations. This is due to the combined model’s noise autocorrelation function displaying a longer correlation length (67 km) in contrast to the RCR model’s (17 km). Consequently, the combined model exhibits a greater reduction in noise variance for height anomaly differences relative to white noise compared to the RCR model.


Introduction
High-resolution (quasi-)geoid models are commonly used on national and continental scales to support GNSS levelling or to replace national vertical reference surfaces based on spirit levelling. However, the quality of these models is often incompletely described, limited to comparisons with geometric height anomalies from GNSS/levelling at a set of height markers, or restricted to information about noise standard deviations based on oversimplified models of noise in the input datasets and simplifications in the process of covariance propagation (e.g., Kearsley 1986;Agren 2004;Voigt et al. 2009;Featherstone et al. 2011;Denker 2013 Brown et al. 2018). A full-noise variance-covariance (VC) matrix of the (quasi-)geoid has not yet been published or analysed, despite its crucial role in providing a comprehensive description of the quality of a (quasi-)geoid model, and its importance as a prerequisite for quality control and further data processing.
In this paper, we aim to address this issue by analysing the full-noise VC matrices of two different quasi-geoid models computed from real data for an area between 2.5 • E and 6.5 • E and 48.5 • N and 56.5 • N, which includes the Netherlands mainland, continental shelf and Wadden islands. The first solution was computed using the traditional removecompute-restore (RCR) technique and is referred here to as the RCR solution. The second solution used a satelliteonly global geopotential model (GGM) with a full-noise VC matrix as an additional dataset alongside the high-resolution datasets used in the computation of the RCR solution. This solution is referred to as the combined solution. We estimated the model parameters of each solution using weighted leastsquares and variance component estimation (VCE) for data weighting. We then computed a full-noise VC matrix of the model parameters using the law of covariance propagation.
Special attention was given to defining the noise VC matrix of all input datasets, including several high-resolution datasets such as airborne gravity disturbances, land and shipboard gravity anomalies, and along-track height anomaly differences from radar altimeter sea surface heights corrected for the instantaneous dynamic topography. For terrestrial gravity anomalies, we extended the functional model by adding a bias parameter per dataset to account for systematic and random errors. The methodology of model parameter estimation for both the RCR and combined solutions was discussed in Klees et al. (2017Klees et al. ( , 2018; the pre-processing of the radar altimeter instantaneous sea surface heights was discussed in Slobbe (2013), the noise model for height anomaly differences from reduced radar altimeter sea surface heights was discussed in Farahani et al. (2017), and the two quasigeoid models considered in this study were presented and analysed in Slobbe et al. (2019).
The remainder of this paper is organised as follows. In Sect. 2, we briefly present the spherical radial basis function models of the residual disturbing potential used in this study and outline the parameter estimation procedure. While the covariance propagation of the RCR solution is straightforward, it is not for the combined solution due to a non-vanishing cross-covariance matrix of the estimated coarse and fine scale model parameters. We derive the expression for the cross-covariance matrix in Sect. 3. In Sect. 4, we summarize the datasets and stochastic models used in this study. We validate the proper scaling of the two noise VC matrices of the residual disturbing potential and derived linear functionals in Sect. 5. Finally, we present the analysis of the full-noise VC matrices of a set of height anomalies in Sect. 6. We conclude this paper by summarising the main findings and discussing the stochastic models of the datasets, which are elaborated on in Sect. 7.

Parameterisation and parameter estimation
The RCR solution uses a single-scale spherical radial basis function (SRBF) model of the residual disturbing potential, Here, "residual disturbing potential" refers to the difference between the Earth's gravitational potential and the sum of normal gravity potential (regularized GOCO05s model com-plete to degree 280) and RTM gravitational potential. The SRBFs {Φ 1 (x, z q,1 ) : q = 1 . . . Q 1 } with centres {z q,1 : q = 1 . . . Q 1 } were chosen to resolve all wavelengths present in the data up to the finest scales resolved by the high-resolution datasets. The model parameters {d q,1 : q = 1 . . . Q 1 } were estimated using weighted least-squares and variance component estimation (VCE) for data weighting. The estimated model iŝ Moreover, a full-noise VC matrix of the estimated model parameters was computed from the noise VC matrices of the input datasets using the law of covariance propagation. The combined solution requires the use of a two-scale SRBF model (cf. Klees et al. 2017), (3) T cs is a low-resolution model representing the coarse scales and T fs adds the finer scales which are not resolved by the coarse-scale model up to the finest scales the high-resolution datasets allow to resolve. For simplicity, we refer to T cs and T fs as the coarse-scale and fine-scale model, respectively. The coarse and fine scale model parameters were estimated sequentially using the methodology in Klees et al. (2017). First, we restored GOCO05s from degree 151 to 280 to the estimated modelT RCR of Eq.
(1) with model parameters {d 1,q : q = 1 . . . Q 1 }, applied a low-pass filter, and synthesized height anomalies on a low-resolution N -point Reuter grid (Reuter 1982) forming the first low-resolution dataset, F i is the height anomaly functional evaluated at point x i , P is the low-pass filter and * is spherical convolution. The associated noise VC matrix was computed by strictly applying the law of covariance propagation. The second low-resolution dataset was computed at the same Reuter grid from the lowpass-filtered GOCO05s model as whereT GGM is the disturbing potential computed from the satellite-only GGM, which now is the difference between the Earth's gravitational potential and the sum of normal gravity potential (regularized GOCO05s complete to degree 150) and RTM gravitational potential. To compute the fullnoise VC matrix of this dataset, we first computed the noise VC matrix of the regularized GOCO05s model from the GOCO05s normal equation matrix N and the inverse of the regularized normal equation matrix, N −1 λ , as N −1 λ N N −1 λ and then applied the law of covariance propagation. Second, the two low-resolution datasets of Eqs. (4) and (5) were used to estimate the parameters {d q,0 : q = 1 . . . Q 0 } of a lowresolution SRBF model of the residual disturbing potential, using weighted least-squares and VCE for data weighting.
Here, Ψ 0 (·, z q,0 ) = P * Φ 0 (·, z q,0 ) is a low-pass-filtered SRBF Φ 0 (·, z q,0 ) with centre z q,0 ; the centres {z q,0 : q = 1 . . . Q 0 } were chosen to resolve all signals present in the low-resolution datasets. The estimated model and parameters are denotedT cs and {d q,0 : q = 1 . . . Q 0 }, respectively. The estimated high-resolution modelT RCR with model parameters {d 1,q : q = 1 . . . Q 1 } was filtered with the high-pass filter H = I − P, where I is an all-pass filter and P is the low-pass filter from Eqs. (4) and (5), to provide an estimate of the fine scales of the residual disturbing potential, Note that the model parameters d 0 and d 1 were always considered as deterministic quantities. Hence, when Tikhonov regularization was applied to the normal matrix, the noise VC matrix of the estimated model parameters was computed as N −1 λ N N −1 λ , where N is the normal matrix and N λ is the regularized normal matrix.
The final estimate of the two-scale SRBF model of the residual disturbing potential follows then from Eq. (3) aŝ withT cs from Eq. (6) andT fs from Eq. (7). The interested reader is referred to Slobbe et al. (2019) for a detailed description of the choice of the SRBFs and the low-pass filter P, the design of the two SRBF networks, the sequential estimation of the two-scale SRBF model, and the results of the VCE.

Covariance propagation
The noise VC matrix for any vector F of linear functionals of the residual disturbing potentialT RCR (Eq. (2)) orT comb (Eq. (8)) can be obtained using the law of covariance propagation applied to FT RCR and FT comb , respectively. While the application of this law is straightforward for FT RCR , it is necessary to take into account that the estimates {d q,0 : q = 1 . . . Q 0 } and {d q,1 : q = 1 . . . Q 1 } are also correlated when applying the law of covariance propagation to FT comb . This correlation arises because the low-resolution dataset {ξ 1 (x i ) : i = 1 . . . N } (Eq. (4)), which is one of the datasets which enters the estimation of the model parameters {d q,0 : q = 1 . . . Q 0 }, also depends on the estimates {d 1,q : q = 1 . . . Q 1 } according to Eqs. (1) and (4).
To derive an expression for the cross-covariance matrix between the two estimatesd 0 andd 1 , we assume that C 0 and C 1 denote their noise VC matrices. We define a set of points {u m : m = 1 . . . M}, e.g., the nodes of an equal-angular grid which covers the domain of the modelT . We are interested in the noise VC matrix of height anomalieŝ The height anomalies of Eq. (9) form a M ×1 vectorξ . Then, using Eqs. (3) and (9), we may writê A is a M × Q 0 matrix, and B 1 , tively. This follows directly from Eqs. (6), (7), and (8). The noise VC matrix of the vectorξ results from the law of covariance propagation applied to Eq. (10): where C 0,1 is the cross-covariance matrix betweend 0 and d 1 and C 1,0 is the cross-covariance matrix betweend 1 and d 0 . With we can write Eq. (11) as The term C cs := AC 0 A T , C fs := (B 1 − B 2 )C 1 (B 1 − B 2 ) T , and C corr := Q + Q T is the contribution of the coarse-scales, the fine-scales, and the covariances between the parameter vectorsd 0 andd 1 , respectively. To compute Q of Eq. (12), we need an expression for the cross-covariance matrix C 0,1 , and for the latter, we need the observation equations for the estimation of the coarse-scale model parameters {d q,0 : q = 1 . . . Q 0 }. They can be written as (cf. Eqs. (4), (5), and (6)) or in matrix-vector notation, The data vector ξ 1 can be written as The elements of the vector ξ 2 are given by Eq. (5). The noise VC matrix of the vector ξ is where C ξ 1 and C ξ 2 is the noise VC matrix of the data vector ξ 1 and ξ 2 , respectively. Then, the weighted least-squares estimated 0 can be written aŝ where is the regularised normal matrix associated with the estimation of the coarse-scale model parameter vector d 0 , R is the regularisation matrix, and λ is the regularization parameter (we assume that regularization was applied; if not, λ = 0). Applying the law of covariance propagation to Eq. (18) and taking into account that Cov(d 1 , ξ 2 ) = 0 provides the expression for the cross-covariance matrix C 0,1 , When using the inversion-free formulation of the weighted least-squares estimator with regularisation as suggested in Klees et al. (2018), Eq. (18) has to be replaced witĥ where λ and λ are regularization parameters (cf. Klees et al. (2018) for the details). To derive the expression for the crosscovariance matrix C 0,1 in this case, we define a matrix which allows writing Eq. (21) aŝ Using the partition where S 1 and S 2 are Q 0 × N matrices, we may writê Then, the cross-covariance matrix C 0,1 follows directly from Eqs. (16) and (25) as Using Eqs. (12) and (26), we can compute the noise VC matrix Cξ of Eq. (13).
We would like to emphasize that all noise VC matrices analysed in this study represent only the effect of random noise and do not include the bias introduced by applying regularisation when estimating the parameter vectors d 0 and d 1 .

Datasets and stochastic models
The RCR model and the combined model were computed for a specific area of interest, which is delineated by the red box in Fig. 1. The figure also displays the different types of local, high-resolution datasets used in the computation, namely 411, 947 along-track height anomaly differences, which were computed from instantaneous sea surface heights from multiple radar altimeter missions corrected for instantaneous dynamic topography from a high-resolution hydrodynamic model (cf. Farahani et al. 2017 for details); 455,335 land surface gravity anomalies; 8205 airborne gravity disturbances; 94,137 shipboard gravity anomalies, and 7179 spatially interpolated gravity anomalies to fill small areas along the Dutch and Belgium coast with no data. The GGM dataset, which is utilized in the computation of the combined model, is based on the regularized GOCO05s spherical harmonic model (Mayer-Gürr et al. 2015). As per the data combination methodology developed in Klees et al. (2017Klees et al. ( , 2018, the GGM dataset comprises a set of height anomalies obtained from a low-pass-filtered GOCO05s model. The low-pass filter is the spectral domain analogue of a Tukey filter with a flat passband from spherical harmonic degrees 2 to 150, a cosine-tapered transition band from degrees 151 to 230, and a stopband above degrees 230. Hence,T cs resolves wavelengths larger than 160 km at the mid-latitudes of the model domain. Wavelength between 110 and 160 km are partially inT cs and partially inT fs due to the transition band of the low-pass filter. Further details about the datasets can be found in Slobbe (2013) 4) and (5), respectively, are located at the nodes of a Reuter grid with grid parameter 240 (Reuter 1982). In Fig. 1, they are displayed as yellow circles marked with a cross. For more information about this choice, we refer to Slobbe et al. (2019). Finally, all datasets were corrected for topographic signals using RTM and harmonic corrections (Forsberg and Tscherning 1981;Heck and Seitz 2007).
In this study, we have employed stochastic models that enhance current practices in several ways. Firstly, we have utilized full-noise VC matrices for the radar altimeter datasets as described in Farahani et al. (2017). Secondly, we have applied a full-noise VC matrix for the GGM dataset, which was derived from the regularized and unregularized GOCO05s normal matrices using the law of covariance propagation, as explained in Slobbe et al. (2019). Thirdly, we have included a bias parameter per land gravity dataset and incorporated the 133 statistically significant bias parameters (at a 5% level of significance) in the final least-squares adjustment. This approach aims to account for systematic and random errors, particularly long-wavelength errors, in the land gravity datasets resulting from issues such as missing atmospheric corrections, height datum inconsistencies, and errors in atmospheric and drift corrections. Fourthly, we estimated variance factors per observation group (i.e., one for land gravity anomalies, marine gravity anomalies, airborne gravity disturbances, interpolated gravity anomalies, and altimeter-based along-track height anomalies, respectively) when computing the RCR model parameters {d q,1 : q = 1 . . . Q 1 }, utilizing MINQUE MCVCE (Kusche 2003). Lastly, we employed MINQUE VCE (Rao 1971) to estimate variance factors for the two coarse scale height anomaly datasets when computing the model parameters {d q,0 : q = 1 . . . Q 0 }.
In all other cases, we used diagonal noise VC matrices. For some of the land gravity datasets noise standard deviations were provided as part of the metadata information, which we adopted as a-priori values. For gravity datasets without metadata information and for the interpolated gravity data, we set the a priori noise standard deviation to a conservative value of 2 mGal. The treated the RTM reductions as deterministic quantities. Future studies may wish to include digital elevation model errors in the stochastic model of reduced datasets.

Validation
Assuming that the noise VC matrices of the input datasets are accurate and parameterization errors are minimal, the law of covariance propagation guarantees that the noise VC matrix of the residual disturbing potential and of any set of height anomalies computed from it such as the matrix Cξ of Eq. (13) is accurate. The largest parameterisation errors are only a few millimeters, and their standard deviation is significantly lower than the standard deviations of the height anomaly noise (as shown in Klees et al. (2017)). However, it is impossible to validate the accuracy of the noise VC matrices of the input datasets. To assess whether the noise VC matrices of the height anomalies evaluated at the GPS height markers are appropriately scaled, we used GPS-levelling data for the Netherlands. The basic idea was to compare the standard deviation of the observed differences between gravimetric and geometric height anomalies with the formal standard deviation obtained from the noise VC matrices of gravimet-ric height anomalies, GPS ellipsoidal heights, and levelled heights.
The histogram in Fig. 2 illustrates the differences between geometric and gravimetric height anomalies at 82 GPS height markers in the Netherlands. The RCR solution exhibits a standard deviation of 0.66 cm, while the combined solution has a standard deviation of 0.95 cm. We compared these values with the values computed from the corresponding noise covariance matrices. During the adjustment of the GPS and levelling network, we found noise standard deviations of approximately 0.5 cm for GPS-ellipsoidal heights and noise standard deviations ranging between 0.2 and 0.5 cm for levelled heights. Consequently, the formal noise standard deviations of geometric height anomalies lie within a range of 0.54 cm to 0.71 cm. The noise standard deviations of gravimetric height anomalies at the GPS height markers fall between 0.07 and 0.44 cm (for the RCR solution) and 0.80 cm to 1.10 cm (for the combined solution) (cf. Fig. 2b, e). Based on these figures, the formal noise standard deviations of the differences between geometric and gravimetric height anomalies are between 0.54 and 0.71 cm (for the RCR solution) and between 0.97 and 1.31 cm (for the combined solution). The observed standard deviation of 0.66 cm for the RCR solution falls well within this range, while the observed standard deviation of 0.95 cm for the combined solution is slightly outside the interval, but not significantly so. These findings suggest that there are no significant scaling errors in the noise VC matrices of height anomalies for either the RCR solution or the combined solution.  Table 1 presents some statistics. The fine scale contribution, (I − P) * T RCR , dominates the height anomaly noise variances of the RCR solution. On average, fine scales account for 73% of the noise variances, while coarse scales, P * T RCR , account for the remaining 27%. Conversely, for the combined solution, the coarse scales,T cs , account for 98% of the noise variances, while to fine scales,T fs , contribute only 2%. Additionally, Fig. 4 and Table 1 indicate that the contribution of the coarse/fine scale correlations to the noise variances of the combined solution is negligible, accounting for only about 0.05% of the noise variances.

Variance-covariance analysis
The height anomaly noise SDs in Fig. 4 and Table 1 illustrate a significant difference between the RCR and combined solutions. Specifically, the noise standard deviations of the coarse scales are much smaller for the RCR solution, averaging just 0.07 cm, as opposed to the combined solution, where they average 0.96 cm. It is worth noting that the two low-resolution height anomaly datasets, ξ 1 (from the low-pass-filtered RCR model, see Eq. (4)) and ξ 2 (from the low-pass-filtered GOCO05s mode, see Eq. (5)), which enter the estimation of the coarse scale model parameter vector d 0 , already exhibit large differences in noise standard deviations. Prior to the application of MINQUE VCE, these values range from 0.03 to 0.13 cm for ξ 1 and from 2.36 to 2.67 cm for ξ 2 (see Fig. 3 for a geographical representation). MINQUE VCE scaling increased the noise VC matrix of dataset ξ 1 by a factor of 509.9, resulting in noise standard deviations now ranging from 0.57 to 2.63 cm. Conversely, the noise VC matrix of dataset ξ 2 was scaled down by a factor of 0.83, leading to noise standard deviations that now range from 2.15 to 2.43 cm.
To better understand this result, it is important to recall that the design matrices of both datasets are identical and square (cf. Klees et al. (2017), Klees et al. (2018)). In other words, if the size of the parameter vector d 0 is denoted by N , then each design matrix has dimensions of N × N . If the noise VC matrices of both datasets were scaled identity matrices, the VCE normal matrix would be singular, making it impossible to estimate variance factors for both datasets (Amiri-Simkooei, 2007). If the noise VC matrices of both datasets were diagonal with little variation along the diagonal, the VCE normal matrix would be ill-conditioned (as considered by Slobbe et al. (2019)). The fact that the VCE normal matrix was well-conditioned and allowed us to estimate two variance factors is due to the fullness of both noise VC matrices.
We also tested other VCE estimators, but the variance factors obtained were very similar to those provided by MINQUE VCE. Therefore, we believe that the most likely explanation for the VCE result is the different structure of the noise autocorrelation matrices of the two datasets, which are shown in Fig 5. The figure indicates that the noise autocorrelation matrix of datavector ξ 2 has more and larger positive and negative noise autocorrelations compared to datavector   ξ 1 , resulting in a stronger weighting of datavector ξ 2 relative to datavector ξ 1 . However, we cannot exclude the possibility that other factors contributed to the unexpected VCE result, such as the extremely high condition number of the a priori noise VC matrix of datavector ξ 2 , which required some regularization to obtain a reasonable solution (as discussed in detail by Klees et al. (2018)). Unfortunately, to our knowledge, there is currently no method available to verify this. Figure 6 displays a plot of the height anomaly noise autocorrelation matrix for the RCR solution and the combined solution, as well as the autocorrelation matrices of height anomaly noise associated with the coarse and fine scales, respectively, at the nodes of an equal-angular grid. It is important to. note that the term "coarse scales" refers to the passband of the low-pass filter P, which was a Tukey filter with a transition band beginning at spherical harmonic degree 150 (equivalent to about 160 km full wavelength at the mid-latitudes of the model domain) and ending at spherical harmonic degree 230 (about 110 km). The half-power-point effective cut-off frequency is at spherical harmonic degree 180, corresponding to a full wavelength of approximately 135 km at the mid-latitudes of the model domain. The pat-tern of the autocorrelation matrix for the RCR solution is dominated by the fine scales, as shown in Fig. 6a and c while the coarse scales dominate the pattern of the autocorrelation matrix for the combined solution, as shown in Fig. 6e and d. This difference can be attributed to the different magnitudes of the noise variances between coarse and fine scales for each solution. Specifically, the ratio of fine scale to coarse scale noise variances is on average 4.4 for the RCR solution, compared to 0.02 for the combined solution. Additionally, Fig. 6 shows that the autocorrelations for the combined solution frequently switch their sign and decrease slowly, while those for the RCR solution do not oscillate strongly and decrease quickly.
A more detailed understanding of the noise autocorrelations can be gained by computing autocorrelation functions. Figure 7a and c show the mean omnidirectional height anomaly noise autocorrelation function (MOACF) over the domain of the modelT comb (indicated by the red box in Fig. 1) and the 1−σ -scatter about the mean for the RCR solution and the combined solution, respectively. The MOACF associated with the combined solution exhibits positive autocorrelations over distances ranging from 0 to 100 km, with quickly Note that the fine scale autocorrelation matrix is the same for the RCR solution and the combined solution by construction decaying oscillating autocorrelations for larger distances. For example, beyond a distance of 220 km, the autocorrelations drop below 0.10. The scatter about the mean reflects heterogeneity and anisotropy of the autocorrelations, never exceeding 15% of the mean.
The MOACF of the RCR solution differs significantly from that of the combined solution. In particular, the autocorrelation length is much shorter, approximately 17 km compared to 67 km for the combined solution, and there are no oscillations beyond the second zero crossing. Additionally, the scatter about the mean is significantly larger for distances below 100 km compared to the combined solution, indicating stronger heterogeneity and/or anisotropy of the autocorrelations (see Fig. 7b).
MOACFs can also be computed for sub-areas. An example is shown in Fig. 8 for the Netherlands and Belgium. A comparison with Fig. 7a and c indicates that they do not differ much from the MOACFs over the domain of the modelT comb (indicated by the red box in Fig. 1). Likewise, the MOACFs for the Netherlands and Belgium are very similar, though the graphs reveal small variations in the location of the first zero-crossing. For the Netherlands, the first zero crossing occurs at 110 km (for both RCR and combined solutions), while for Belgium it is located at 87 km (for the RCR solution) and 100 km (for the combined solution). Moreover, slightly larger negative autocorrelations between the first two zero-crossings are observed for Belgium compared to the Netherlands (for both RCR and combined solutions).
The autocorrelations of the noise have implications for GNSS-levelling, as positive values indicate lower noise standard deviations for height anomaly differences compared to white noise, and vice versa. Figure 9 illustrates the reduction in percent of the standard deviation of height anomaly differences above white noise, as a function of the distance between two points in the Netherlands and Belgium, respectively. The plot reveals significant reductions in noise standard deviations for the combined solution, and smaller reductions for the RCR solution, depending on the distance. Table 2 presents the reduction percentages for selected distances. For example, the combined solution for the Netherlands provides a reduction of at least 86% for distances not exceeding 10 km. The reduction is still greater than 53% for distances not exceeding 40 km. Between the first two zero-crossings, located at about 100 km and 230 km, respectively, there is a slight degradation due to negative noise autocorrelations, which does not exceed 15%. Similar results are obtained for  Belgium. The RCR solution shows much smaller reduction percentages. For instance, over the Netherlands, the reduction is only 39% for a distance of 10 km and 22% for a distance of 40 km. For Belgium, the reductions are somewhat larger, but still small compared to those of the combined solution. Evidently, the use of GOCO05s as one of the noisy datasets, in combination with a stronger weight assigned to this dataset compared to the high-resolution datasets by VCE, is responsible for the more favorable noise autocorrelation function of the combined solution compared to the RCR solution. Instead of expressing the reduction of noise in height anomaly differences as a percentage, we can also compute the noise standard deviations directly as a function of the distance between two points. Figure 10 provides an example of this approach. Over the Netherlands, the combined solution shows very low noise standard deviations of only a few millimeters over distances of a few tens of kilometers. How-ever, they increase with distance, reaching over 1 cm at about 70 km. For the RCR solution, the noise standard deviations are even smaller, not exceeding 0.2 cm for distances shorter than 80 km. Similar result are obtained for Belgium.

Discussion and conclusions
We conducted an analysis of the full-noise VC matrices for two regional quasi-geoid models that cover Belgium and the Netherlands, including the mainland, continental shelf, and Wadden islands. The official quasi-geoid model of the Netherlands was computed using the RCR technique from land, marine, and airborne gravity and radar altimeter data. The second model used the same datasets as the RCR model and a satellite-only GGM as an additional noisy dataset.
For the RCR solution, we found that the height anomaly noise standard deviations range between 0.03 cm and 0.85 cm Fig. 10 Mean noise standard deviation (SD) of height anomaly differences as function of the distance between two points (blue) compared to white noise (red) for the Netherlands (a, c) and Belgium (b, d). RCR solution (a, b) and combined solution (c, d) over the domain of the modelT comb . Fine scales (i.e., wavelength shorter than about 110 km at mid-latitudes of about 52 • ) accounted for about 73% of the noise standard deviations. In contrast, the height anomaly noise standard deviations for the combined solution ranged between 0.53 cm and 1.58 cm. Here, coarse scales dominated, accounting for approximately 98% of the noise standard deviations.
Our analysis showed that the noise VC matrix's structure for height anomalies is favorable for GNSS-levelling, especially for the combined solution. The noise standard deviations of height anomaly differences are significantly smaller than expected from the height anomaly noise standard deviations due to positive autocovariances. Height anomaly differences' precision improves as the distance between two points decreases. For the combined solution, we obtained standard deviations of less than 0.2 cm for distances below 10 km and less than 0.6 cm for distances below 40 km over the Netherlands mainland. We also found that the favorable positive autocovariances result from using the GOCO05s satellite-only GGM as one of the noisy datasets when computing the combined solution.
Additionally, we observed that the RCR solution provided a mean omnidirectional height anomaly noise autocorrelation function with a correlation length that was almost four times shorter. The benefit of positive noise autocorrelations when computing height anomaly differences is less pronounced. Nonetheless, the noise standard deviations of height anomaly differences were still smaller than those of the combined solution.
The quality of the full-noise VC matrix of a quasi-geoid model relies on various factors, including the quality of the stochastic models, which cannot be validated. In our study, we made significant efforts to ensure reasonable VC matrices for all involved datasets, especially for the GGM dataset and the radar altimeter datasets. Additionally, we accounted for long-wavelength systematic and random errors in the land gravity datasets through bias parameter estimation. We employed VCE for weighting among land, marine, airborne, and interpolated gravity, and radar altimeter datasets, and when combining them with the GGM dataset. Furthermore, we strictly followed the law of covariance propagation while estimating the model parameters.
The most critical simplification in our stochastic model is neglecting existing noise correlations in land and marine gravity datasets. This omission could significantly affect the results of data combination, and provide over-optimistic noise VC matrices for both the RCR and the combined models. It is highly probable that noise in shipboard gravity data is correlated due to errors in the drift and Eötvös corrections, errors in the connection with land gravity stations, and the use of cross-over adjustment. Although there are no studies on noise correlations in marine point gravity anomalies, some studies have been published on noise in marine mean gravity anomalies. For instance, Monka et al. (1979) analyzed 1 • ×1 • mean free-air gravity anomalies in the North Sea region and reported error covariances of 3 mGal 2 at distances of 1000 to 2000 km. Weber and Wenzel (1983) analysed 6 × 10 mean free-air sea gravity anomalies in the North Atlantic, North Sea, and Mediterranean Sea, and reported error covariances of more than 5 mGal 2 at distances up to 100 km and still 3 mGal 2 at distances up to 1, 000 km. Noise correlations are also likely to exist in land gravity data. Nonetheless, even if noise correlation models exist in the land and marine gravity datasets, the numerical complexity of utilizing this information in the least-squares adjustment is considerable, given the size of these datasets.
Another way to enhance the stochastic model is by taking into account any errors in the digital elevation model (DEM) and the mass densities utilized in the computation of the RTM reduction. However, as mentioned earlier, it remains uncertain if such information can be incorporated into the least-squares adjustment process, as it entails noise correlations across all datasets, leading to even greater numerical complexity.
Using another satellite-only GGM which is provided with a full-noise VC matrix will change the results for two reasons. Firstly, its noise VC matrix will differ from the one used in this study. Secondly, the parameters of the low-pass filter which was used when combining the land, marine, and airborne gravity and radar altimeter datasets with the GGM dataset may change. However, it is very likely that the favourable structure of the noise VC matrix of height anomalies of the combined solution compared to a RCR solution will persist for any satellite-only GGM which uses data of the satellite gravity missions CHAMP, GRACE, GOCE, and GRACE-FO.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no relevant financial or nonfinancial interests to disclose. All authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript. The authors have no financial or proprietary interests in any material discussed in this article. The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.