A methodology for least-squares local quasi-geoid modelling using a noisy satellite-only gravity field model

The paper is about a methodology to combine a noisy satellite-only global gravity field model (GGM) with other noisy datasets to estimate a local quasi-geoid model using weighted least-squares techniques. In this way, we attempt to improve the quality of the estimated quasi-geoid model and to complement it with a full noise covariance matrix for quality control and further data processing. The methodology goes beyond the classical remove–compute–restore approach, which does not account for the noise in the satellite-only GGM. We suggest and analyse three different approaches of data combination. Two of them are based on a local single-scale spherical radial basis function (SRBF) model of the disturbing potential, and one is based on a two-scale SRBF model. Using numerical experiments, we show that a single-scale SRBF model does not fully exploit the information in the satellite-only GGM. We explain this by a lack of flexibility of a single-scale SRBF model to deal with datasets of significantly different bandwidths. The two-scale SRBF model performs well in this respect, provided that the model coefficients representing the two scales are estimated separately. The corresponding methodology is developed in this paper. Using the statistics of the least-squares residuals and the statistics of the errors in the estimated two-scale quasi-geoid model, we demonstrate that the developed methodology provides a two-scale quasi-geoid model, which exploits the information in all datasets.


Introduction
In this paper, we investigate the combination of a noisy satellite-only global gravity model (GGM) with noisy highresolution datasets (e.g. terrestrial gravity anomalies) to estimate a local quasi-geoid model using weighted leastsquares techniques. By considering the satellite-only GGM as one of the noisy datasets, we expect to improve the quality of the estimated local quasi-geoid model. By exploiting existing information about the noise variances and covariances in combination with weighted least-squares techniques, we aim at making a step forward towards a comprehensive description of the quality of the estimated quasi-geoid model in terms of a full noise covariance matrix for quality control and further data processing.
The problem is timely. The quality and spatial resolution of the most recent satellite-only GGMs, which are mainly based on data of the Gravity Recovery and Climate Experiment (GRACE) and Gravity field and steady-state Ocean Circulation Explorer (GOCE) satellite missions, have improved dramatically compared to pre-mission models. Moreover, the GGM's spherical harmonic coefficients are now complemented with a full noise covariance matrix. For models such as GOCO05s (Mayer-Gürr et al. 2015), the quality of the noise covariance matrix benefits from (i) a post-fit residual analysis, which provides more realistic models of the data noise (e.g. Farahani et al. 2013), and (ii) numerically efficient algorithms to propagate the full data noise covariance matrices into the noise covariance matrix of the estimated spherical harmonic coefficients (e.g. Farahani et al. 2013). Some progress has also been made regarding high-resolution datasets frequently used in local quasi-geoid modelling. An example is Farahani et al. (2017) who derived coloured noise models for radar altimeterbased along-track quasi-geoid height differences, or Slobbe (2013) who successfully accounted for long-wavelength errors in terrestrial gravity anomalies (e.g. Heck 1990) by augmenting the functional model with additional parameters.
Until now, a GGM with full noise covariance matrix has not been used as one of the noisy datasets in the computation of a local quasi-geoid model. The standard approach is remove-compute-restore approach, i.e. the GGM is used to facilitate a local approach to quasi-geoid modelling by removing most of the energy in the data at the long wavelengths. From a theoretical point of view, the combination of a local set of terrestrial gravity anomalies with a noisy GGM has been considered in a number of publications since the 1980s. They all use a modification of Stokes' formula and go back to the early work of Sjöberg and Wenzel (Sjöberg 1980(Sjöberg , 1981Wenzel 1981). At that time, only error degree variances of the GGM were used. The modification of Stokes' kernel was formulated as a global optimization problem, in which either the variance or the mean square error of the quasi-geoid heights was minimized. Since then, the stochastic spectral combination methods were applied routinely when computing local quasi-geoid models. An example is the European quasi-geoid EGG08 (Denker et al. 2008;Denker 2013), which uses the spectral combination method of Wenzel (1981). Sjöberg (2005) was the first who derived the formalism of a local least-squares modification of Stokes' formula, which uses, among others, the full noise covariance matrix of the GGM to compute the weights per spherical harmonic degree. In (Sjöberg 2011), the restriction to weights per spherical harmonic degree was given up. To our knowledge, this method has not been investigated yet in detail. Some numerical aspects were studied in Ågren (2004) and Ellmann (2004), in particular numerical instabilities when estimating the spectral weights, which naturally arise in local applications.
Here we follow another approach to local quasi-geoid modelling, which uses least-squares techniques to estimate the parameters of a local spherical radial basis function (SRBF) model of the disturbing potential from the available datasets. Least-squares techniques have many advantages. In particular, they allow us to include noise covariance matrices, improve them using variance component estimation, provide variance-covariance information about the estimated parameters and linear functionals of them, and allow the use of statistical hypothesis testing to test the validity of the mathematical model. Least-squares local quasi-geoid modelling using SRBFs has been intensively studied by various authors, see Klees et al. (2008) for a literature overview until 2007. Since then, a number of studies investigated various aspects of the use of SRBFs in local quasi-geoid modelling. This comprises the choice of the type of SRBFs (e.g. Tenzer and Klees 2008;Bentel et al. 2013a, b), SRBF network design and numerical optimization (e.g. Wittwer 2009), regularization issues (e.g. Naeimi 2013), and the optimization of the location of the SRBFs (e.g. Lin et al. 2014). Some aspects related to the combination of data with different bandwidths have been discussed in Panet et al. (2011);Naeimi (2013); Bentel and Schmidt (2016); Lieb et al. (2016);Lieb (2017). However, they do not cover numerical studies about the combination of a GGM with full noise covariance matrix with high-resolution noisy datasets.
An alternative to the use of a GGM as one of the noisy datasets in local quasi-geoid modelling is to complement the high-resolution local datasets with the original satellite data at altitude, e.g. satellite gravity gradients from the GOCE mission, low-low satellite-to-satellite tracking data from the GRACE mission, and high-low satellite-to-satellite tracking data from GRACE, GOCE, and other low Earth orbiters. Examples in the context of local quasi-geoid modelling with SRBFs are Lieb et al. (2015) and Lieb (2017). The major drawback of this approach is the complexity of the functional model for the low-low satellite-to-satellite tracking data and the huge amount of data. This may be the reason why numerical studies published so far (e.g. Lieb 2017) do not use high-low satellite-to-satellite tracking data, limit to a subset of the available GOCE data (e.g. the second radial derivative of the gravitational potential), or use GRACE-based alongtrack gravity potential differences as pseudo-observations instead of the original K-band ranging data, often in combination with simplified noise models. Overall, this approach does not offer a significant advantage compared to the use of a GGM that is based on the same data and, therefore, is not pursued in this study.
The study addresses to main research questions: (i) What is a suitable functional model for the satellite-only GGM? (ii) How to combine the satellite-only GGM with high-resolution datasets to obtain a quasi-geoid model that optimally exploits the information content in all datasets?
The remainder of this paper is organized as follows. First, we introduce local single-scale and two-scale SRBF models of the disturbing potential and suggest three functional models to be used in a least-squares estimation of the quasi-geoid model from the satellite-only GGM and the high-resolution datasets. Following this, we describe the setup of the numerical experiments, which were designed to investigate the performance of the SRBF models and the functional models. Thereafter, we present and discuss the results of the numerical experiments. We conclude by emphasizing the main findings and identifying topics for future research.

Parameterization and functional models 2.1 Local parameterization of the disturbing potential
A prerequisite for local gravity field modelling is that the involved datasets do not contain much long-wavelength energy, where "long" relates to the size of the data area. One possibility to achieve this is by reducing all datasets for the contribution of a global model of the disturbing potential. Then, the disturbing potential to be parameterized in local modelling is a residual quantity with little, though nonzero energy at the long wavelengths. Moreover, there is an upper limit of the highest attainable spatial resolution, which depends on the data distribution and signal-to-noise ratio. This allows us to consider a local model of the disturbing potential which is band-limited to some maximum degree, say L max . From now on, we call this residual and band-limited disturbing potential simply "disturbing potential", denoted T .
In this study, we use SRBFs to model T over an area of interest. Basically, two models will be used. The first one is a single-scale model, i.e. where is a SRBF located at z i , {c i } are the SRBF coefficients, which are to be estimated from the data using least-squares techniques, Q l is the reproducing kernel of the space of harmonic functions of degree l, φ l is the Legendre coefficient of degree l,x andẑ i are points on the unit sphere, and σ R is the surface of a sphere of radius R. The model of Eq. (1) is referred to as a single-scale model. Alternatively, we may use a multi-scale model involving several sets of SRBFs representing different scales, i.e.
where j indicates the scale, {c j· } are the SRBF coefficients at scale j, and Ψ j (·, z ji ) is a SRBF of scale j centred at the point z ji . In the context of a multi-scale analysis, the SRBF Ψ j may be defined as (e.g. Lieb 2017) Frequently, the relation l j = 2 j − 1 is used to relate the scale index j to the maximum spherical harmonic degree l j , which is resolved at scale j, though other choices are possible.

Functional models
In the framework of this study, we assume that there are basically two datasets, i.e. a low-resolution dataset and a high-resolution dataset. The low-resolution dataset d 1 is synthesized from the spherical harmonic coefficients of the GGM as where {ĉ nm } are the spherical harmonic coefficients of the GGM, {c (ref) nm } are the spherical harmonic coefficients of the reference GGM, and H nm is a solid spherical harmonic of degree n. The low-resolution dataset is band-limited to a degree L 1 ≤ L GGM , where L GGM is the maximum degree of the GGM. We assume that the high-resolution dataset {d 2 (x 2k ) : k = 1 . . . K 2 } allows the resolution of wavelengths up to a maximum degree L 2 ≤ L max , where L 2 depends on the point density and the signal-to-noise ratio. Defining a kernel a spherical convolution of T with δ L as and linear functionals F 1 and F 2 of the disturbing potential T , we may relate the datasets d 1 and d 2 to the disturbing potential T as where E{·} denotes mathematical expectation. We will investigate three functional models to estimate a local quasi-geoid model by least squares from the low-resolution dataset d 1 and the high-resolution dataset d 2 .
Functional model no. 1 uses the single-scale model of the disturbing potential, Eq. (1), and reads: The coefficients {c i } are estimated simultaneously from the two noisy datasets using weighted least-squares techniques.
The weight matrix of each dataset is the inverse of the noise cofactor matrix.
Functional model no. 2 uses the single-scale model of the disturbing potential, Eq. (1), and reads The kernel P of Eq. (12) is defined as The Legendre coefficients {h n : n = 1, 2, . . .} are equal to 1 for degrees n ≤ p 1 , taper off between degrees p 1 < n < p 2 , and are zero for all degrees n ≥ p 2 . An example is a cosine taper, This taper will be used in the numerical experiments of Sect. 3. The coefficients {c i } are estimated simultaneously from the two noisy datasets using weighted least-squares techniques. The weight matrix of each dataset is proportional to the inverse of the noise covariance matrix. The noise covariance matrix of P * d 1 is computed from the full noise covariance matrix of d 1 using the law of covariance propagation.
The difference between the functional models no. 1 and no. 2 is in the functional model of the low-resolution dataset. Functional model no. 2 uses a tapered SRBF, whereas functional model no. 1 uses a truncated SRBF. Moreover, functional model no. 2 applies the same taper to the dataset, whereas functional model no. 1 uses the original dataset.
Functional model no. 3 uses a two-scale model of the disturbing potential, i.e. Eq. (3) with J = 2: The first term on the right-hand side is a low-resolution model of T comprising degrees from 0 to L 1 , i.e. its resolution is identical to the resolution of dataset d 1 . The second term on the right-hand side complements the lowresolution model to the maximum resolution L 2 of dataset d 2 . In the context of a multi-resolution analysis, it represents a detail space comprising wavelengths from degrees L 1 + 1 to L 2 . The basis functions Ψ 1 and Ψ 2 of Eq. (16) are defined as with Φ(x, z) of Eq.
(2). Inserting the last two equations into Eq. (16), the two-scale model of the disturbing potential T is written as with Φ of Eq.
(2). The coefficients {c 1i } and {c 2i } are estimated in two steps. First, we use the functional model and estimate the coefficients {c 2i } using weighted least squares. Then, we define a new dataset where {ĉ 2i } denotes the least-squares estimate of {c 2i }. The resolution of the dataset d 3 is identical to the resolution of the dataset P * d 1 . In that sense, d 3 and P * d 1 are spectrally consistent. Then, we use the functional model and compute an estimate {ĉ 1i } of the coefficients {c 1i }, using weighted least-squares techniques. The noise covariance matrix of dataset d 3 is computed from the noise covariance matrix of the estimated coefficients {ĉ 2i } using the law of covariance propagation. It is a full matrix. The least-squares estimate of the disturbing potential is then given by Eq. (19), with {c 1i } and {c 2i } replaced by the estimates {ĉ 1i } and {ĉ 2i }, respectively.

Remarks
1. The motivation of using the functional model of Eq. (10) is the following. Dataset d 1 and its full noise covariance matrix are band-limited to a degree L 1 ≤ L GGM . Therefore, the right-hand side of the functional model must also be band-limited to the same degree. To achieve this, we consider the signal 1·Φ, with Φ of Eq.
(2). We expand this signal on the sphere σ R in spherical harmonics and truncate the expansion at degree L 1 . The result is identical to δ L 1 * Φ. If the right-hand side of the functional model would not be band-limited to degree L 1 , the leastsquares estimate of the coefficients {c 1i } would be biased towards zero for the wavelengths above degree L 1 . This is due to the fact that a band-limited noise covariance matrix is equivalent to zero noise and noise correlations for degrees above L 1 . 2. The motivation to use the functional model of Eq. (12) is the result of numerical experiments which are described in Sect. 3 and discussed in Sect. 4. There, we will show that |d  (2017), which in our notation is Moreover, Lieb (2017) suggests to estimate the coefficients {c 1i } and {c 2i } simultaneously using weighted least-squares techniques. Some preliminary experiments with this model and Φ set equal to the Abel-Poisson kernel (Freeden et al. 1998) point to a sub-optimal quality of the estimated quasi-geoid model at the resolution of the dataset d 1 , which is likely caused by the simultaneous estimation of the two sets of coefficients {c 1i } and {c 2i }. However, additional numerical experiments are necessary to support these preliminary results. They are out of the scope of this study.

Numerical experiments
The parameterizations and functional models of Sect. 2 will be analysed using numerical experiments. Though from a practical point of view, working with real data may be desired, we decide to use a state-of-the-art combined GGM and a satellite-only GGM to generate the (noise-free) highresolution and low-resolution datasets, respectively. The main motivation for us to prefer GGMs to real datasets is that some problems and limitations of the functional models of Sect. 2.2 would be masked by deficiencies in real datasets, e.g. unmodelled signal and noise and data gaps. This would make a proper interpretation of the results impossible. Generating the exact data from GGMs and adding noise which is consistent with the corresponding noise covariance matrix provides a complete error control and facilitates a proper interpretation of the results. The datasets are generated from the GGMs using a spherical harmonic synthesis and thereafter reduced for the contribution of a long-wavelength GGM, which serves as the reference model. Here we use EIGEN-6C4 (Förste et al. 2014) to generate the (noise-free) high-resolution dataset and the regularized version of GOCO05s (Mayer-Gürr et al. 2015) to generate the (noise-free) low-resolution dataset. The latter is also used as the reference model, though up to a smaller maximum degree. The low-resolution dataset consists of a set of height anomalies. This is a logic choice as the target quantity is a quasi-geoid model. The noise covariance matrix of the low-resolution dataset is obtained from the full noise covariance matrix of the spherical harmonic coefficients of the unregularized version of GOCO05s by applying the law of covariance propagation. A logic choice for the  (20): Poisson wavelets of order 3, depth = 30 km, Fibonacci grid with a mean node distance of 25 km; max absolute parameterization error = 4.5 mm, SD of parameterization error = 0.9 mm (over area of interest) Ψ 1 of Eq. (22): Poisson wavelets of order 3, depth = 60 km, Fibonacci grid with a mean node distance of 60 km; max absolute parameterization error = 0.01 mm (over area of interest) The area of interest is 49 • -63 • N and 6 • W-10 • E in all set-ups The parameterization area is identical to the data area high-resolution dataset would be gravity anomalies. Here, we use gravity disturbances for simplicity reasons. Noise in gravity disturbances is zero-mean white Gaussian, i.e. the noise covariance matrix of the high-resolution dataset is a scaled unit matrix. More details about the datasets are provided in Table 1. The datasets are generated at the Earth's surface and cover an area, which is referred to as "the data area". The Earth's surface is represented by the digital elevation model Euro-DEM v1.0 (Hovenbitzer 2008) with 2 grid width. In areas where this model is not available, we use SRTM version 2.1 (Farr et al. 2007) with 3 grid width. For the remaining areas, we use ASTER GDEM v2 (Tachikawa et al. 2011) with 1 grid width.
The SRBF of Eq.
(2) is a Poisson wavelet of order 3 (Holschneider and Iglewska-Nowak 2007). Different from the Shannon kernel, which is frequently used in local quasigeoid modelling, the Legendre spectrum of a Poisson wavelet relative to a sphere of radius R has a peak at degree 3R R−|z| , where z < R is the location of the Poisson wavelet (cf. Fig 1). The Legendre spectrum of a Poisson wavelet may give the impression that a single-scale Poisson wavelet model is not able to accurately represent a quasi-geoid with a resolution one typically encounters in practice. Therefore, Chambodut et al. (2005) suggest to use Poisson wavelets of different scales to guarantee that the space of spherical harmonic complete to a degree L 2 is sufficiently well covered. However, Slobbe (2013) successfully used a single-scale Poisson wavelet model to compute a quasi-geoid model for the Netherlands mainland, continental shelf, and Wadden Islands with an accuracy of about 1.5 cm standard deviation using real data. The only prerequisite is that the energy in the data at the lowest and highest frequencies is reduced by using a reference GGM and a digital terrain model, respectively. The Poisson wavelets are located at a constant depth below the Earth's surface and cover the data area. Their horizontal positions correspond to the points of a Fibonacci grid (Gonzalez 2010). Whenever a new set of Poisson wavelets is chosen in the numerical experiments, we have to determine the optimal depth and the optimal mean distance between the Poisson wavelets. This is done using noise-free datasets generated from the corresponding GGMs on grids dense enough to preserve the information content in the GGM. We define a set of candidate depths and candidate mean distances, and estimate the model coefficients by least squares using the corresponding dataset (i.e. gravity disturbances when looking for a highresolution model and height anomalies when looking for a low-resolution model). The depth and mean distance that provide the model with the smallest RMS difference to a height anomaly control dataset are selected. The fit of this model to the control dataset is referred to as the "parameterization error". Note that the parameterization error is always defined in terms of height anomalies, no matter whether the dataset comprises gravity disturbances or height anomalies. Other control datasets are generated to assess the quality of the estimated quasi-geoid models in Sect. 4. They always comprise height anomalies on grids different from the data grids and are computed using a spherical harmonic synthesis of the GGMs from which the noise-free datasets were generated. When computing quasi-geoid models using weighted least-squares techniques, we calculate the normal equations explicitly and apply Tikhonov regularization (Tikhonov 1963) with a unit regularization matrix. The regularization parameter is fixed using the method in (Wittwer 2009). The normal equations are solved using a parallelized QRdecomposition with column pivoting. This solver is preferred to a Cholesky decomposition due to its much better stability for ill-conditioned linear systems at the benefit of a smaller bias in the least-squares estimate due to a smaller regularization parameter. Table 1 summarizes the set-up of the numerical experiments, which will be used in Sect. 4.

Functional model no. 1
We use experimental set-up no. 1 of Table 1. Table 2 shows some statistics of the least-squares residuals for the estimated quasi-geoid model. The standard deviation (SD) of the residuals is 9.86 cm for the low-resolution dataset and 2.02 mGal for the high-resolution dataset. The latter corresponds to the standard deviation of the noise in the high-resolution dataset. The former, however, is much larger than the noise. From this we conclude that the estimated quasi-geoid model fits the high-resolution dataset within noise, but gives a poor fit to the low-resolution dataset. Table 3 shows the statistics of the errors in the estimated quasi-geoid model. They are computed over the area of interest.
The errors range from −10.68 to 12.76 cm; the error SD is 3.46 cm. After applying a low-pass filter at the cut-off degree 200, the estimated quasi-geoid model error SD is 7.62 and 7.68 cm depending on what low-resolution signal is taken as the reference.
To get more insight into the reason why the estimated quasi-geoid model does not fit the low-resolution dataset within noise, we repeat the experiment with noise-free data. The error SD of the estimated quasi-geoid model reduces from 3.46 to 0.09 cm. This is identical to the SD of the parameterization error. Hence, when using noise-free data, the esti- The statistics are computed over the area of interest, which is bounded by 49 • -63 • N and 6 • W-10 • E, on an equal-angular grid of width 10 48 (1st row) and 23 29 (2nd and 3rd rows) T 1 = GOCO05s minus GOCO05s complete to degree 150; T 2 = EIGEN-6C4 minus GOCO05s complete to degree 150. F 1 = height anomaly functional mated quasi-geoid model perfectly fits the high-resolution dataset. This does not apply, however, to the low-pass-filtered quasi-geoid model; the error SD is 7.12 and 7.46 cm, respectively, i.e. comparable to the results using noisy data. Hence, the poor fit of the estimated quasi-geoid model to the lowresolution dataset cannot be explained by the noise in this dataset. Additional numerical experiments (not shown here) reveal that the fit to the low-resolution dataset can only be improved by further increasing the size of the data area. The 5 • extension beyond the area of interest in all directions as used here is already a challenge in real quasi-geoid modelling as access to data of neighbouring countries is not guaranteed. Moreover, we found that the fit to the low-resolution dataset improves slowly when enlarging the data area. From this we conclude that the poor fit of the estimated quasi-geoid model to the low-resolution dataset is caused by the hard truncation of the Poisson wavelets. This introduces strong spatialdomain oscillations, which are cut off at the border of the data area when computing the elements of the design matrix. This introduces errors in the functional model, which exceed by far the noise in the low-resolution dataset, as shown in Table 3.

Functional model no. 2
We use experimental set-up no. 2 of Table 1. The kernel P, which according to Eq. (12) is used in the functional model of the low-resolution dataset, is chosen according to Eqs. (14) and (15). The cosine taper parameters are set equal to p 1 = 150 and p 2 = 230. Hence, the filtered low-resolution dataset P * d 1 of Eq. (12) is band-limited to degree L 1 = 230. The choice of p 1 and p 2 is a trade-off between loss of information in the low-resolution dataset by filtering (i.e. nonzero d 1 − P * d 1 ), and a reduction in the area under the side lobes of the cosine taper, which cause oscillations of the filtered Poisson wavelets extending beyond the data area. Moreover, the difference p 2 − p 1 determines how fast the oscillations roll off. The difference p 2 − p 1 = 80 has been fixed after some numerical experiments. Note that the maximum degree of the reference GGM (which is L ref = 150 in our experiments) and the maximum degree of the noisy The statistics are computed over the area bounded by 44.5 • -67.5 • N and 11.5 • W-14.5 • E # Of points refers to the whole dataset GGM (which is L GGM = 280 for GOCO05s) impose lower and upper bounds, respectively, on the choice of p 1 and p 2 , i.e. p 1 ≥ 150 and p 2 ≤ 280. The maximum possible value of p 2 is equal to the maximum degree of the GGM, L GGM . The GGM used in this study is GOCO05s. However, for GOCO05s the cumulative height anomaly commission error increases exponentially with increasing degree. It is 1.5 cm at degree 200, but already 3.6 cm at degree 230, and 6.8 cm at degree 250. The height anomaly signal and noise degree variances intersect at degree 257. Hence, when assuming that the noise standard deviation in the high-resolution dataset does not exceed 1-2 mGal (which applies to good terrestrial gravity anomaly datasets), it does not make sense to use a low-resolution dataset complete to the maximum possible degree of L GGM = 280. Another reason in favour of a choice L 1 < L GGM is the fact that the condition number of the noise covariance matrix of the low-resolution dataset, which is propagated from the full noise covariance matrix of the spherical harmonic coefficients of the GGM, increases with increasing L 1 . This makes the computation of the least-squares estimator numerically challenging. Table 4 shows the statistics of the least-squares residuals for the estimated quasi-geoid model. The SD is 1.94 mGal for the high-resolution dataset. This is close to the SD of the superimposed zero-mean white Gaussian noise of 2.0 mGal. From this we conclude that the model fit is within noise. The situation is different for the low-resolution dataset. The SD of the residuals is 3.42 cm. This is a factor of 2 larger than the average SD of the data noise. Obviously, for some reason, the estimated quasi-geoid model does not fit the low-resolution dataset as one may expect. Table 5 shows the statistics of the errors in the estimated quasi-geoid model. They are computed over the area of interest.
The error SD is 3.01 cm. The error SD of the cosinetapered quasi-geoid model is 1.89 and 1.76 cm, respectively, depending on what control dataset is used. We may compare this with a quasi-geoid model which is estimated using only the high-resolution dataset. The corresponding error SDs are 3.01, 1.58, and 1.71 cm, respectively. Hence, adding the low-resolution dataset does not improve the accuracy of the estimated quasi-geoid model at low frequencies (i.e. below The statistics are computed over the area of interest, which is bounded by 49 • -63 • N and 6 • W-10 • E, on an equal-angular grid of width 10 48 (1st row) and 23 29 (2nd and 3rd rows) T 1 = GOCO05s minus GOCO05s complete to degree 150; T 2 = EIGEN-6C4 minus GOCO05s complete to degree 150. F 1 = height anomaly functional degree 230). This is unexpected, because the low-resolution dataset has a high quality and should improve the estimated quasi-geoid model at the low frequencies. This, and the large least-squares residuals of the dataset P * d 1 , which by far exceed the noise, implies that for some reason, the information content in the low-resolution dataset is not fully exploited in the combined least-squares adjustment.
To help understand this result, we run two additional experiments. First of all, we compute the error of the functional model of the (noise-free) low-resolution dataset, Eq. (12) We found that it does not exceed 0.01 cm over the area of interest. In a second experiment, we compute a low-resolution quasi-geoid model using the noisy low-resolution dataset and the functional model of Eq. (12). The SD of the residuals is 0.13 cm. This is significantly smaller than for the solution which uses both datasets (SD = 3.42 cm, cf. 1st row in Table 4). Moreover, the error SD of the low-resolution quasi-geoid model is 1.58 cm when evaluated over the area of interest. This is also smaller than the error we obtain when using both datasets (SD = 3.89 cm, cf. 1st row in Table 5).
Our interpretation of the results of these experiments is that a single-scale model is not able to fit two datasets of significantly different bandwidths. Consequently, the weighted least-squares principle forces the solution to match the high-resolution dataset (because it comprises many more observations than the low-resolution dataset) at the price of a larger mismatch to the low-resolution dataset.
To support this interpretation, we choose experimental setup no. 3 of Table 1, which is similar to experimental set-up no. 2, but involves a high-resolution dataset and a low-resolution dataset with a much larger bandwidth difference of 335% compared to 117% of experimental set-up no. 2. The SD of the least-squares residuals of the low-resolution dataset increases from 3.42 cm (cf. Table 4, 2nd row) to 5.13 cm, whereas the SD of the least-squares residuals of the highresolution dataset does not change. Hence, when increasing the bandwidth difference between the high-resolution and the low-resolution datasets, the fit of the model to the low-  (20), the dataset d 3 of Eq. (21) would suffer from edge effects. Therefore, the high-resolution dataset d 2 must be available over an area, which is larger than the data area of the dataset d 3 . The additional extension must be chosen to reduce the edge effects in dataset d 3 below the noise level. In the experimental set-up no. 4, we use an additional extension by 5 • in all directions. Test computations reveal that this choice causes edge effects in dataset d 3 , which are negligible compared to the noise. We expect that the additional extension can be chosen much smaller. To find the minimum extension may be the subject of another study. Table 6 shows the statistics of the least-squares residuals for the model of Eqs. (20) and (22), respectively. The fit of the high-resolution dataset d 2 to the model of Eq. (20) has a SD of 1.74 mGal. This is close to the superimposed noise of SD = 2 mGal; a similar fit has also been observed when using the functional model no. 2 (cf. Table 4). However, compared to the functional model no. 2, the fit of the dataset P * d 1 to the model has improved dramatically: from 3.42 cm (1st row of Table 4) to 1.47 cm (2nd row of Table 6). A value of 1.47 cm is consistent with the noise in the dataset P * d 1 . The fit of the dataset d 3 to the model is 2.69 cm, i.e. dataset P * d 1 has a larger contribution to the model compared to dataset d 3 .
(Both datasets are evenly large.) This is also consistent with the expectations based on an analysis of the noise covariance matrices of the two datasets (not shown here).
The two-scale model appears to have a much higher quality than the single-scale model of Sect. 4.2. This follows from The low-pass filter is a cosine taper with parameters p 1 = 150 and p 2 = 230. The statistics are computed over the area of interest, which is bounded by 49 • -63 • N and 6 • W-10 • E, on an equal-angular grid of width 10 48 (1st row) and 23 29 (2nd and 3rd rows) T 1 = GOCO05s minus GOCO05s complete to degree 150; T 2 = EIGEN-6C4 minus GOCO05s complete to degree 150. F 1 = height anomaly functional the statistics of the differences at the control datasets, which are shown in Table 7. For instance, the fit of the two-scale model to the control dataset F 1 (δ 500 * T 2 ) improves from SD = 3.01 cm (2nd row of Table 5) to SD = 1.97 cm (2nd row of Table 7). The fit to the low-resolution control data improves dramatically, too: from 1.89 cm (F 1 (P * T 1 ), 1st row of Table 5) and 1.76 cm (F 1 (P * T 2 ), 3rd row of Table 5) to 0.86 and 0.77 cm, respectively. From this we conclude that the two-scale model in combination with the functional model of Eqs. (20), (22) performs better at all wavelengths than any of the two single-scale models. The improvement is a factor of 2.2 for the wavelengths common to the highand the low-resolution dataset, and a factor of 1.5 for the wavelength not resolved by the low-resolution dataset. The former is due to the fact that the suggested approach which uses the two-scale model fully exploits the higher accuracy of the low-resolution dataset, which is not the case if any of the single-scale models is used.

Summary and conclusions
In this study, we investigated different approaches to estimate a local SRBF model of the disturbing potential using weighted least squares from a high-resolution dataset and a low-resolution dataset. In practice, the low-resolution dataset represents a satellite-only spherical harmonic model of the global gravity field equipped with a full noise covariance matrix. Considering the latter as one of the noisy datasets in local quasi-geoid modelling is considered as a significant improvement to the traditional remove-computerestore approach. It improves the quality of the estimated quasi-geoid model and paves the way to a complete quality description of the estimated quasi-geoid model in terms of a full noise covariance matrix. Two approaches investigated in this study use a singlescale SRBF model, but differ in the functional model for the low-resolution dataset. The third one uses a two-scale SRBF wavelet model and estimates the coefficients per scale independently of each other.
We showed that the functional model of the low-resolution dataset has to be chosen with care. A hard truncation of the SRBFs at the maximum degree of the low-resolution dataset is the right choice in global quasi-geoid modelling, but provides a wrong functional model in local quasi-geoid modelling. This is in line with the results in (Slobbe et al. 2012). Applying a taper to both the low-resolution dataset and the SRBF model solves this problem.
We also showed that a single-scale SRBF model cannot deal with datasets of different bandwidths. The estimated quasi-geoid model is biased towards the high-resolution dataset at the cost of a poor fit to the low-resolution dataset. The latter appeared to be much worse than the noise in this dataset suggested, which indicates that the information content of the low-resolution dataset is not fully exploited.
We suggested the use of a two-scale SRBF model in combination with a sequential estimation of the scale-dependent coefficients. The latter differs from what has been suggested in the literature in the context of a multi-scale analysis. In this way, we ensure that the two datasets are weighted in line with their accuracy, the information content in the low-resolution and high-resolution datasets is fully exploited, and the misfit of the estimated quasi-geoid model is consistent with the noise in the datasets.
A challenge of the suggested approach in applications involving real datasets is the additional extension of the data area for the high-resolution dataset. In this study, a safe choice has been made to make edge effects insignificant. In applications involving real datasets, access to high-quality terrestrial gravity anomaly datasets of neighbouring countries is not guaranteed. How much the data area needs to be extended and whether data with reduced accuracy can be used in the additional area without introducing distortions in the estimated quasi-geoid model has to be investigated.
It would be interesting to compare the two-scale approach suggested in this study with a multi-scale approach, which estimates the coefficients at the two scales simultaneously as suggested in Lieb (2017). Some preliminary experiments (not shown here) indicate that such a multi-scale approach provides a sub-optimal low-resolution quasi-geoid model compared to a sequential estimation as suggested here. Whether this may be corrected for by a further optimization of the multi-scale approach, for example, by introducing constraints between the model coefficients associated with different scales may be the subject of a future study.