Integration of airborne gravimetry data filtering into residual least-squares collocation: example from the 1 cm geoid experiment

Low-pass filters are commonly used for the processing of airborne gravity observations. In this paper, for the first time, we include the resulting correlations consistently in the functional and stochastic model of residual least-squares collocation. We demonstrate the necessity of removing high-frequency noise from airborne gravity observations, and derive corresponding parameters for a Gaussian low-pass filter. Thereby, we intend an optimal combination of terrestrial and airborne gravity observations in the mountainous area of Colorado. We validate the combination in the frame of our participation in ‘the 1 cm geoid experiment’. This regional geoid modeling inter-comparison exercise allows the calculation of a reference solution, which is defined as the mean value of 13 independent height anomaly results in this area. Our result performs among the best and with 7.5 mm shows the lowest standard deviation to the reference. From internal validation we furthermore conclude that the input from airborne and terrestrial gravity observations is consistent in large parts of the target area, but not necessarily in the highly mountainous areas. Therefore, the relative weighting between these two data sets turns out to be a main driver for the final result, and is an important factor in explaining the remaining differences between various height anomaly results in this experiment.


Introduction
In this paper we adapt the residual least-squares collocation (RLSC, Willberg et al. 2019) so that correlations from low-pass-filtered airborne gravity observations are handled consistently. Simultaneously, we present our final result from an International Association of Geodesy (IAG) joint working group (JWG) which is designed to support the realization of the International Height Reference System (IHRS, Ihde et al. 2017). Within this JWG 2.2.2, called 'the 1 cm geoid experiment', different participating groups calculate height anomaly, geoid height and potential values by using identical terrestrial and airborne gravity observations. The main objective of this JWG is to increase the compatibility between different methods for regional geoid determination B Martin Willberg martin.willberg@tum.de 1 Institute of Astronomical and Physical Geodesy, Technical University of Munich, Arcisstrasse 21, 80333 Munich, Germany by analyzing and quantifying differences between the results. Furthermore, the experiment should help to define common standards for the IHRS realization and verify the quality of the submitted results. Additional information about the purpose and benefit of the JWG is published in a summary paper , this issue).
The study area of 'the 1 cm geoid experiment' lies at the southern end of the Rocky Mountains in the United States of America; mainly in the states of Colorado and New Mexico. It was selected for its geoid slope validation survey from 2017 (GSVS17) where positions, gravity and deflections of the vertical are measured with very high accuracy at a line of 223 benchmarks along U.S. Highway 160. However, these measurements or their processing are not yet published, so they will function as reference values for the JWG results only in the future. Furthermore, the region has good coverage in terms of terrestrial and airborne gravity measurements, but is also intended to be a very challenging region for regional geoid determination, as it includes highly mountainous areas. Another major challenge within the JWG, and one main topic of this paper, is the optimal inclusion of airborne grav-ity observations from the 'Gravity for the Redefinition of the American Vertical Datum' project (GRAV-D, GRAV-D Team 2018b). Within this project, the complete area of the United States is covered with equally distributed airborne observations in order to define a new gravity-based vertical datum, preferably with an accuracy of 1-2 cm (GRAV-D Team 2017).
For our calculation we include RLSC (Willberg et al. 2019) as a modified version of the standard least-squares collocation (LSC) by Moritz (1980). By including a removecompute-restore (RCR) approach with a high-resolution global geopotential model (GGM) and a topographic gravity model, the input of RLSC consists only of residuals. As GGM, we include the XGM2018 model, an internal successor of XGM2016 (Pail et al. 2018. Both models provide a full variance-covariance information from regional varying weighting, which enables improved error-modeling in RLSC. However, XGM2018 provides a more realistic stochastic model for the GGM reduction (see discussion on this issue in Willberg et al. 2019). Additionally, accuracy estimation in RLSC benefits from the fact that individual covariance matrices contain stochastic information of all involved components. Detailed justification for the application of RLSC and a comparison between RLSC and standard LSC is published in Willberg et al. (2019). This contribution adapts RLSC to airborne gravimetry for the first time.
A crucial task thereby is to separate the target gravity signal from the observation noise which is normally implemented by a low-pass filter (e.g., Childers et al. 1999). While the variety of different filters for this purpose is huge, their purpose is very similar: to reduce the highfrequency noise from the airborne observations (details in Sect. 2). The application of a low-pass filter is necessary for platform-stabilized gravimeters (Childers et al. 1999;Olesen 2003) as well as strapdown gravimeters (Wei and Schwarz 1998;Becker 2016;Li 2011). Detailed descriptions of the two different measurement systems and the predominant error sources are given in Schwarz and Wei (1995). In general, both gravimeter types are used within the GRAV-D project, but in the target area (5th block mountain south, MS05), we only have airborne measurements from the platform-stabilized gravimeters Micro-g LaCoste TAGS S-137 (turn-key airborne gravimetry system) and TAGS S-211. Detailed documentation about the instrumentation and the flight lines from block MS05 is provided in GRAV-D Team (2018a).
The inclusion of a low-pass filter for the processing of airborne gravity data is the standard procedure. However, a low-pass filter will inevitably result in a significant correlation of the airborne observations along the track of the aircraft, which is usually not considered in airborne processing. As an example, Forsberg et al. (2000Forsberg et al. ( , 2014 and Hwang et al. (2007) assume low-pass filtered airborne observations to be uncorrelated in LSC. Within this paper, for the first time, we derive an approach which includes correlations resulting from a low-pass filter consistently in the functional and stochastic model of RLSC. As a result of the modification, we have a more consistent error modeling and filter dependent covariance matrices. Taking these correlations into account seems of even higher significance when the airborne observations are combined with other measurements (e.g., terrestrial gravity observations).
We see another advantage of our presented approach in the fact that the combination of terrestrial and airborne gravity measurements is included in a direct, one-step LSC, which also contains field transformation and the downward continuation of airborne measurements. For the combination of different data sets, the one-step calculation allows the full exploitation of the much higher airborne resolution in alongtrack direction. For other approaches, which either calculate a regular grid before the downward continuation (Forsberg et al. 2000) or include a spectral method for the analysis of airborne measurements at flight height (Smith et al. 2013;Jiang and Wang 2016), the difference between along-track and across-track accuracy gets lost.
In summary, the main innovations of this paper include (1) a novel approach to include a low-pass filter to regional gravity field modeling in general, and RLSC in particular; (2) application of GRAV-D data with recommendations for its use; (3) the first application of RLSC to real gravity observations, as the theory paper on RLSC (Willberg et al. 2019) uses a synthetic simulation environment; (4) demonstration of the RLSC contribution to 'the 1 cm geoid experiment' and comparison of corresponding results. This paper is structured as follows. In Sect. 2 we motivate the application of a low-pass filter for airborne observations. Next, in Sect. 3 a suitable filter is derived and included into the RLSC formalism, resulting in a rigorous formulation for the combination of filtered (airborne) and unfiltered (terrestrial) observations. The data sets and our calculation procedure from 'the 1 cm geoid experiment' are explained in Sect. 4, and its results are analyzed in Sect. 5. Additionally, we include comparisons of our results along GSVS17 in relation to the JWG mean value (Sect. 6). Lastly, we draw conclusions and give an outlook in Sect. 7.

Reasons for low-pass filter in airborne gravimetry
For optimal comparisons within 'the 1 cm geoid experiment' the input observations used by all contributing groups should be identical. Therefore, the airborne observations of block MS05 are provided as gravity measurements within the JWG. In contrast to the original GRAV-D airborne data available at the homepage of the National Geodetic Survey (2019, NGS), the provided data is already corrected for the individual line bias by comparison with already existing models and crosspoint validation. Furthermore, the data is provided with the uniform frequency of 1 Hz while the gravimeter Micro-g LaCoste TAGS S-211 originally records data at 20 Hz (GRAV-D Team 2018a). As is well known, airborne gravity observations include in general higher noise levels than stationary terrestrial data. This results from the dynamics of the aircraft. Corrections have to be applied for Eötvos accelerations (Harlan 1968), the vertical acceleration (Zhong et al. 2015), and the offlevel error of the instrument (LaCoste 1967;Swain 1996). The exact corrections that are applied to the observations are described in GRAV-D Team (2017), and a description of the underlying software with the use of the same notation can be found in Zhong et al. (2015). It is commonly assumed that the occurring noise shows up mainly in the short wavelengths, while the gravity signal dominates the longer wavelengths (e.g., Childers et al. 1999). Accordingly, the inclusion of a low-pass filter in along-track direction of airborne gravimetry is common practice.
Airborne observations from the GRAV-D project are already low-pass filtered with a time domain Gaussian filter (GRAV-D Team 2017, Chapter 2.2). However, its purpose was to preserve the amplitude of the gravity signal, and it leaves significant short-wavelength noise in the data (GRAV-D Team 2018a). Accordingly, the GRAV-D Team (2018a, Chapter 4.1) recommends using a second low-pass filter for the airborne observations in order to remove the short-wavelength noise. In this study, we assume that the first filter actually leaves the gravity signal untouched (as it was intended), because it allows us to neglect the first filter in the theoretical derivations (Sect. 3). As a consequence, in the following we call the GRAV-D observations provided within 'the 1 cm geoid experiment' as 'original airborne observations', although they already include the first low-pass filter, the line bias correction and the resampling to uniform 1 Hz frequency. Accordingly, in further usage of the expression "filter", we are referring only to the second one. The GRAV-D Team (2018a, Chapter 4.1) also emphasizes that the low-pass filter should be included prior to a downward continuation, as it would otherwise amplify the noise. Since LSC implicitly includes a downward continuation, we apply a suitable low-pass filter before the LSC method (details in Sect. 3).
In general, we consider the characterization of the highfrequency noise as a sophisticated problem, because its amplitude and frequency distribution can change due to various reasons like wind conditions or flight velocity. An analysis of the most dominant error sources in airborne gravimetry and its assessment can be found in Schwarz and Wei (1995). However, we need to specify filter criteria in order to separate the noise and gravity signal. We demonstrate this in the following example: some of the flight lines in MS05 are actually divided into two segments (GRAV-D Team 2018a). Accordingly, there is a total of seven flight segments, where the end of one flight line overlaps with the beginning of another. One of these flight segments has overlapping measurements over approximately 100 km which is significantly more than all other combinations. In Fig. 1, we present the overlapping segments of these two flight lines (FL) 103 (red) and 203 (blue) with its original gravity disturbances (bright solid lines) and the filtered observations (dark solid lines), respectively. Additionally, the dotted lines result from an approximation by available spherical harmonic (SH) coefficient models (details in Sect. 4). Therefore, in our case we include a combination of XGM2018 and the topographic gravity model dV_ELL_Earth2014, which is calculated according to Rexer et al. (2016), but continued to SH degree 5480 by including 1 arc-min topography information from Hirt and Rexer (2015). Note that the presented flight lines result from the same flight plan, which is why their horizontal offset is negligible. Their gravitational difference from the vertical offset is small and can be approximated by the difference between the two modeled observations (dotted lines).
In this section of the flight lines we see major anomalies due to the fast-changing topography, and the gravity disturbances change with about 100 mGal over a distance of only 50 km. The two close-up analyses in Fig. 1 display that there is basically no correlation between the high-frequency signals of FL103 and FL203, which indicates dominating noise in these frequencies. The modeled airborne observations contain gravity signals to SH degree 5480. Nevertheless, they appear very smooth in comparison to the original gravity observations, which is assumed to result from high-frequency noise in the observations as well. By including a Gaussian low-pass filter we reduce high frequencies in the observations, so that the modeled observations and the filtered observations approximately contain the same spectral signal content (Fig. 1). Details about the Gaussian filter are explained at a later stage in Sects. 3.1 and 3.3. Apart from the reduction of high-frequency noise and the recommendations by the GRAV-D team, we consider a low-pass filter necessary for the following reasons: (1) The original airborne observations are available with a frequency of 1 Hz, all in all resulting in 283,716 observations. Because there is a very high correlation between consecutive observations, we do not consider it reasonable to include all observations directly in the processing. However, decreasing the amount of observations is only possible if a low-pass filter is applied first, and otherwise we will create aliasing errors (details in Sect. 3.3).
(2) Furthermore, decreasing the amount of airborne observations is necessary in order to prevent RLSC from containing numerically singular covariance matrices. These would inevitably arise when airborne observations are included with a 1 Hz sampling frequency. (3) In the current implementation our RLSC approach for airborne observations is limited to the maximum SH degree 5480, as this is the highest available resolution for a topographic gravitational potential model (in our case: dV_ELL_Earth2014). Consequently, airborne gravity signals above this SH degree would result in a model error. It is commonly assumed that airborne observations do not include a significant gravity signal above SH degree 5480 (approx. 4 km in spatial resolution) due to signal attenuation with altitude. However, we have to consider the following: a few of the flight lines are located at 5200 m above the ellipsoid while the Rocky Mountains beneath the corresponding lines reach heights of more than 4000 m. Accordingly, the flight tracks can be quite close to the signal generating masses. At the same time, there might be significant errors in the topographic gravity model due to problems of Shuttle Radar Topography Mission (SRTM) in areas with fast-changing topography. Lastly, the limited spatial resolution of topographic elevation models in highly mountainous areas is another reason that might result in errors. Only by including a low-pass filter can we guarantee that these problems do not arise in (R)LSC.

Methodology
In this section, we show the attenuation of different lowpass filters in the frequency domain by using the recursion formulas by Jekeli (1981). Next, we derive a Gaussian lowpass filter that is suitable for airborne gravity data, include it consistently in the functional and stochastic model of RLSC, and explain its effect. Thereby, we also present a numerical example for the final covariance matrix. Lastly, we formulate the combination of filtered airborne and unfiltered terrestrial gravity observations.

A Gaussian low-pass filter
A normalized Gaussian weighting function w G can be written as a function of the Euclidean distance d whereby the standard deviation of the Gaussian filter σ defines the smoothing effect. The weighting function w G (d) includes the normalization factor α so that the overall signal content is not amplified or weakened by the filter. This normalization factor α will be derived later on (cf. Eq. 9). Jekeli (1981) derives the weighting function w G of a Gaussian filter on the sphere in dependence of its spherical distance ψ where b gives the smoothing effect which is defined by the radius of the Earth R Note that Eqs. 1 and 2 can be derived from each other using simple trigonometry. Furthermore, Jekeli (1981) gives recursion formulas which are normalized on the unit sphere. They enable the formulation of the attenuation β, which results from a Gaussian filter per SH degree n in the frequency domain We see that with decreasing β, the attenuation of the Gaussian filter increases when either the SH degree or the standard deviation σ becomes larger. Strictly speaking, the recursion formulas are only valid for points on the sphere, which is why Jekeli (1981) includes a spherical approximation for points on the ellipsoid. By analogy, we include the same approximation since we apply the filter to airborne observations with a relatively constant flight height above the ellipsoid. In the following, the standard deviation σ is replaced by the more intuitive half width at half maximum (HWHM) The HWHM of a Gaussian low-pass filter describes the distance after which the weighting function w G decreases to half its maximum. Figure 2 presents the attenuation β per SH degree n (Eq. 4) that is evoked by different Gaussian filters with a HWHM between 1 and 5 km. While the HWHM of 5 km (green) basically erases all signals at SH degree 5000, the HWHM of 1 km (blue) attenuates only 20% of the signal at this SH degree.

Including the Gaussian filter into the functional model
The Gaussian weighting function defined in Eq. 1 can be formulated to describe the filter process in form of a functional model A G . It calculates filtered observations l G from Note that number and position of the original observations l and the filtered observations l G do not have to be identical. If there are, for example, less point positions in l G than there are in l, the functional model A G will result in a reduction of the sampling frequency. This process is called downsampling in the following. In order to keep the signal content consistent before and after the filter process, A G is normalized so that the summation of one row in A G equals one N j=1 a i j = 1.
Therein, a i j are the matrix elements of the functional model A G and N is the number of observations l in one flight line. Consistently,å i j are the elements of the functional model before the normalization withẘ G (d i j ) being the related Gaussian weighting function for the distance d i j . The normalized Gaussian weighting function can be formulated by adding the normalization factor α i to Eq. 1 We see that the elements of A G after the normalization are defined by the Gaussian weighting function w G (d i j ). Note that the high-frequency noise of airborne observations occurs only in along-track direction, so that we use a 1-D Gaussian filter along individual flight lines. Accordingly, our approach differs from the 2-D filters in Jekeli (1981) or Huang et al. (2008), who apply and normalize their filters on the unit sphere or surface element, respectively. In the functional model A G , every row refers to a filtered observation whereas the columns refer to the original observations. Accordingly, one row in A G defines which original observations are used to calculate the corresponding filtered value. One column in A G describes to what extent a specific original observation will be included in different filtered values (Eq. 6). Visualizing the values of either one row or one column would look like a bell-shaped curve since the functional model A G is calculated from a (Gaussian) low-pass filter. Its shape is thereby defined by the combination of the filter radius in the Gaussian weighting function and the flight velocity at recording time. Note that in general the functional model is not symmetric since it is normalized only in one dimension.

Spectral limitation of airborne gravity data
In Sect. 2 we found that the very high frequencies in the airborne observations should be removed before the downward continuation. Furthermore, the total number of observations has to be reduced before LSC. We regard the functional model of a Gaussian filter A G , as a suitable tool to put these two aspects into practice.
The airborne gravity measurements are discrete observations of the continuous gravity field. According to the Nyquist-Shannon sampling theorem, a continuous signal with the maximum frequency f max can only be reconstructed from an equidistant sampling with a frequency higher than 2 f max . Accordingly, simply reducing the sampling frequency of very high-frequency observations generally results in aliasing effects. In practice, a common approach to treat this issue is the use of a low-pass filter (Childers et al. 1999). It reduces the maximum frequency f max and therefore also the new sampling frequency that is needed to fulfill the Nyquist-Shannon sampling theorem. In the case that the equidistant sampling after the filter process is done with the frequency 2 f s , all frequencies below f s could be restored while the frequencies above f s are irretrievably lost.
We want to exploit these aspects in order to separate signal and noise in the airborne observations. For this purpose, we assume that the airborne observations mainly include noise above a specific SH degree μ and follow our three-step procedure in order to remove this noise by means of the Gaussian functional model A G . Additionally, Fig. 3 visualizes the three steps with their relationship in terms of SH degree and attenuation factor β (introduced in Eq. 4 and Fig. 2).

Fig. 3
Relationship between the SH degree and the attenuation β of the Gaussian filter with μ set to SH degree 5400. The red 'areas' are removed by the Gaussian filter. The two dark areas are assumed to mainly consist of noise, while the two brighter colored areas are dominated by signal (1) A Gaussian low-pass filter removes the two red parts in the frequency domain (Fig. 3). Thereby, it filters most of the noise, but also parts of the signal from the observations. The Gaussian filter radius is selected as a compromise between the preservation of the signal and the removal of the noise. In Fig. 3 we set it exemplary to HWHM = 3 km, as this keeps only 10% of the noise at SH degree μ.
(2) Reducing the sampling frequency of the observations to SH degree μ irretrievably removes the previously filtered data above this SH degree ( Fig. 3: dark red area). By doing this, the dark blue area is discarded as well, however it is still present in the filtered observations, thus causing a small aliasing effect. For the implementation we combine step 1 and 2, so that the functional model A G includes the low-pass filter and the downsampling process (refer to Sect. 3.2). (3) Lastly, the previously filtered signal representing SH degrees below μ ( Fig. 3: light red area) can be restored from the filtered observations l G by means of a highpass filter. This step works analogous to a deconvolution by multiplication with the inverse Gaussian functional model A S . It results in observations l S which are spectrally limited to SH degree μ Here, the Gaussian functional model A S is a square matrix that describes the point positions of l G in both dimensions and is normalized according to Eq. 10. Note that the Gaussian functional models A G and A S are derived in the same way. They differ only in the corresponding point positions, and the normalization factor, as it is related to the input positions.
Applying these steps removes high-frequency noise and results in spectrally limited observations l S which could be inserted into LSC. In this case, the spectral limitation from step one to three could be interpreted as an airborne preprocessing. However, in the next section we include these three steps consistently within RLSC instead, thereby considering the correlations from the Gaussian filter between consecutive observations as well. Linking the described steps to the RLSC procedure, we insert the filtered observations l G (step 2) as input and apply a high-pass filter (step 3) implicitly within RLSC.

Adapting RLSC to include the Gaussian filter model
The RLSC formula for calculating an output s from the input observations l is given by Willberg et al. (2019) where C ll is the covariance matrix of the observations, Cˆlˆl the covariance matrix of an unbiased reduction modell, and Cˆsˆl the cross-covariance matrix between input and output.
In general, the hat-operator marks quantities which have been derived from external reference models (GGM, topographic gravity model). The positions and functionals of covariance matrices are defined by its subscript. In Eq. 12 a remove-compute-restore (RCR) approach analogous to Willberg et al. (2019) is already included, so that we reduce the observations l with unbiased model observationsl before RLSC. The reduction model in the output functionalŝ is restored accordingly after RLSC. Note that the covariance function of the reduction model Cˆlˆl is contained consistently in Eq. 12. Following Eq. 6 we apply the Gaussian functional model A G to l and, due to linearity of this operator, correspondingly to both l andl The corresponding covariance functions of l G andl G are calculated from simple error propagation. From Eq. 6, we obtain the covariance matrix of filtered observations C G and accordingly the covariance function of the model obser- Similarly, we obtain for the cross-covariance matrix C Ĝ considering that the output modelŝ is not filtered, thus multiplying A G only once. From Eqs. 14 to 17, we conclude that the RLSC method with filtered observations is identical to the standard formula for RLSC (Eq. 12) except that the covariance matrices C G , the observations l G and the model observationsl G refer to filtered values instead of the original ones. Note that the output s should contain an unfiltered signal, therefore we do not use superscripts for s andŝ.
In the case that airborne gravity observations are included in Eq. 18, their spectral limitation (Sect. 3.3, step 1 and 2) is already included in the filtered values l G andl G . The deconvolution (step 3) happens within the multiplication of the inverted covariance matrices C G ll and C Ĝ ll .

A noise covariance function for the observations
In LSC, the covariance function of the observation noise is often included as a diagonal matrix (Moritz 1980;Arabelos and Tscherning 2009) under the assumption that different observations are not correlated among each other. Obviously, this assumption does not hold for low-pass-filtered observations as they contain a strong correlation among consecutive measurement epochs. In this section, we derive a covariance function for the filtered observations C G ll in Eq. 15. The covariance information of the unfiltered observations C ll is in general unknown, which is one of the main reasons for applying the Gaussian low-pass filter in the first place. Thus, instead of propagating C ll to C G ll we can define the variances σ G l • σ G l (from elementwise multiplication •) on the main diagonal of C G ll . Thereby, we assume to know the accuracy of the observations after the filter process. In this case C ll is introduced as identity matrix I, and the row vector m is included as factor to set the variances of C G ll to σ G l • σ G l . We can write Fig. 4 Visualization of the covariance function C G ll from FL103 whereby a presents the covariance function for all 4099 observations, while b is a zoom to the first 300 observations whereby Eq. 20 uses elementwise division and diag(X) defines the main diagonal of a square matrix X. By including the identity matrix in Eq. 19 we assume an uncorrelated error for the airborne observations l before the Gaussian filter process which, in general, is not the case. However, the simplification is justified in our case as the correlation of this error is assumed to be small in comparison to the correlation resulting from the Gaussian filter.
Exemplary, Fig. 4a shows the covariance function C G ll (Eq. 19) of a single flight line (FL103) without downsampling. The main diagonal of C G ll is defined by the variances σ G l • σ G l of the filtered observations, which are equally set to 1 mGal 2 . The same information with more detail is presented in Fig. 4b which zooms to the first 300 observations. We can see a band matrix with very strong correlations among consecutive observations instead of the usual diagonal matrix (e.g., Forsberg et al. 2000Forsberg et al. , 2014Hwang et al. 2007). Thus, Fig. 4b verifies that the reduction of the sampling frequency is reasonable in this case due to the very strong correlations among consecutive observations. Furthermore, this is necessary as in general the covariance matrix C G ll without downsampling tends to be numerically singular for the same reason. The downsampling procedure we introduced in Sect. 3.2 would generally retain the structure of a band matrix in C G ll , but the covariance would drop to zero much faster.

Combining filtered and unfiltered observations
In the Colorado experiment, we have a combination of terrestrial and airborne gravity observations, and all reasons we listed for applying a low-pass filter to the airborne observations (Sect. 2) are not valid for the terrestrial ones. Therefore, we present a formulation that applies the standard RLSC from Eq. 12 to the terrestrial observations, while the airborne observations are processed according to Eq. 18 with a Gaussian filter. In order to realize a consistent formulation, we keep naming the original observations l, but now con-sider them as a combination of terrestrial l ter and airborne l air observations Similarly, we handle the reduction modell l = l ter l air .
The combination of filtered and unfiltered measurements requires updating the functional model A G . We now name it A g , as it contains the Gaussian filter only for the airborne part A G is calculated according to Eq. 10 and I is an identity matrix with a size that equals the number of terrestrial observations. We multiply the new functional model A g to the observations l, the model observationsl and the covariance matrices Cˆsˆl and Cˆlˆl. Accordingly, we obtain the RLSC formula for combined input quantities analogous to Eq. 18 In this form, the covariance matrix C g ll defines the input observation accuracy. It is set up as a combination of an unfiltered terrestrial covariance matrix C ter,ter and a filtered airborne covariance matrix C G air,air since the observation noise is not correlated between airborne and terrestrial measurements. The terrestrial part simplifies to (σ 2 ter I) in case the observations are uncorrelated and of equal accuracy. Thereby, σ 2 ter is the variance of the terrestrial observations. The variance σ G l • σ G l of the airborne observations is defined as main diagonal of C G ll (Eq. 19). In the case they are of equal accuracy, σ G l • σ G l can be replaced by the scalar values σ 2 air as simplified variance of the airborne observations. Analogous to Eq. 24, we derive the stochastic model for a combination of filtered and unfiltered observations. The error covariance matrix E ss of the output s is based on the standard formula by Willberg et al. (2019) Therein, Cˆsˆs is the error covariance function for the model observations, which stays unfiltered as it refers only to the output functional.

Data and calculation
From Sect. 2, we conclude that the airborne gravity observations have to be low-pass-filtered in along-track direction. In Sect. 3 we derive a suitable concept to do so. Originally, there were 283,716 airborne gravity observations from the GRAV-D block MS05. They were corrected for their individual line biases and provided with a uniform sampling rate of 1 Hz. We chose a Gaussian low-pass filter with a HWHM of 3 km for its good compromise between removing noise and keeping gravity signal. Fig. 3 demonstrates that this HWHM removes approximately 90% of the noise above SH degree 5400. Additionally, Fig. 1 shows that the resulting frequencies in the filtered airborne observations are similar to the modeled observations.
A consequence of the low-pass filter is the very high correlation between consecutive observations (Fig. 4), which allows a significant reduction of the sampling frequency. However, this reduction is another benefit-risk assessment: a strong downsampling increases the numerical efficiency, but could also result in a loss of signal information. We exemplary show this in Table 1 for five different sampling frequencies and three different filter lengths (given in HWHM). The table gives the mean error in mGal when all 283,716 filtered observations are reproduced from the downsampled observations with a simple spline interpolation. We interpret this value as a quality criterion for the signal loss due to the low-pass filter. Considering that the maximum error of a specific observation can be much higher than the average value, Table 1 Overview of the benefit-risk assessment between HWHM and sampling frequency Colored values are classified by their average signal loss [mGal], which indicates the ability to restore the filtered but non-downsampled signal again. Additionally, combinations are marked which result in covariance matrices with an especially high condition numbers ('bad κ') or even numerical 'singularity' we want to keep the average error well below 0.1 mGal. Furthermore, we mark combinations which are not feasible due to a very high condition number κ for the covariance matrix C G ll respectively with 'bad κ' or numerical 'singularity'. While for a higher sampling frequency or a larger filter length the condition number κ of the covariance matrix C G ll increases, a smaller sampling frequency or shorter filter length increases the signal loss.
We conclude, that first, downsampling is needed in any case to prevent a numerical singular covariance matrix C G ll . Secondly, the sampling frequency of 1/32 Hz is the best compromise for a HWHM of 3 km in the target area. Note that this sampling frequency is a main driver for the final filter characteristics. Together with the flight velocity, it defines the SH degree μ above which we assume mainly noise in the airborne observations (Fig. 3). By reducing the sampling frequency to 1/32 Hz, we result in a total number of 8976 filtered airborne observations. Thereby, two adjacent observations have an average distance of 3.4 km in along-track direction. The actual distance or spatial sampling, however, depends on the flight velocity of the aircraft and varies significantly in the target area. In general, the aircraft have been slowest in the mountainous areas of MS05 so that the shortest point distances are in those areas where it is most beneficial. Otherwise, we could have included a different downsampling method which we tested as well: it provides a constant spatial sampling and therefore results in a varying sampling frequency. However, in this case the filter method would model space-correlated noise in the airborne observations. Since the airborne gravity observation noise is mainly time-correlated instead, our method with constant sampling frequency stays preferable.
The 59  Input residuals of RLSC: reduced gravity disturbances for a the terrestrial observations and b the airborne observations. The observations are removed by the XGM2018 model and the topographic gravity potential. The extent of the output grid is marked by a red box and GSVS17 by a black line are densely measured regions in the target area, but also some observation gaps. The quality of the terrestrial observations is not well-defined as there is almost no metadata available. However, some additional information about the terrestrial data set can be found in Wang et al. ((2020), this issue). The topographic heights of the output area grid (red box in Fig. 5) are visualized in Fig. 6.
In RLSC we include an a priori assumption for the accuracy of the observations (Eq. 25 or Willberg et al. 2019) which is basically unknown for both the terrestrial and the airborne observations. GRAV-D Team (2018a) contains an error analysis for the airborne measurements, but we do not think it can give us reliable accuracies for the following reasons: (1) the derivation of the error is based on observations from different flight heights. Therefore, all values are continued to a mean height by standard free-air correction (GRAV-D Team 2017) which introduces assumptions and errors. (2) The quality assessment in GRAV-D Team (2018a) still contains the high-frequency noise which was reduced in Sect. 3. Lastly, (3) the crossover statistics cannot describe the overall accuracy accurately (GRAV-D Team 2018a), as they only assess measurement errors in the data lines (in our case: east-west), not the crossover lines (north-south).
For the filtered airborne observations we assume a uniform standard deviation of σ air = 1 mGal, which seems reasonable according to Wei (1995), Childers et al. (1999), Novák et al. (2003) and Lu et al. (2017). While the airborne observations are well-distributed, terrestrial observations in general benefit from measuring the full gravity frequency spectrum on the Earth's surface. Therefore, we prefer a RLSC combination where both airborne and terrestrial observations have a similar influence on the overall result. In our case, this can be achieved by a uniform standard deviation of σ ter = 3 mGal for the terrestrial observations, since their number is significantly higher, and they do not include corre- A remove-compute-restore procedure is essential for RLSC and explained in detail in Willberg et al. (2019). In it, we include XGM2018 as high-resolution GGM to its maximum SH degree 760. The topographic gravitational potential model dV_ELL_Earth2014, or the ERTM2160 gravity model (Hirt et al. 2014) are added, respectively, in the frequencies above. Following aspects should be noted regarding the remove-compute-restore procedure: (1) Although, the maximum SH degree of XGM2018 is 760, its actual spectral resolution is defined by the maximum degree in spheroidal harmonics, which is 719. The same effect is visible in dV_ELL_Earth2014 where the maximum SH degree is 5480, but 5400 in spheroidal harmonics, respectively. For simplicity, we stick to the models actual spectral resolution, and as of now the corresponding SH degrees refer to spheroidal harmonics.
(2) Furthermore, the reduction models for terrestrial and airborne data are not identical in our case. ERTM2160 is given in an approximately 250 m resolution grid format, therefore providing the very high-resolution information which we need for the reduction of terrestrial data on the Earth's surface. However, dV_ELL_Earth2014 is calculated from SH synthesis which can be used for arbitrary point positions (e.g., airborne observations). In this case, the different reduction steps are justified as dV_ELL_Earth2014 and ERTM2160 are both calculated from the SRTM v4.1 topography model by Jarvis et al. (2008) in the target area. Accordingly, they contain identical signal information in the respective SH degrees. (3) Lastly, dV_ELL_Earth2014 and ERTM2160 are based on simplified density assumptions and cannot reflect lateral density variations. Deviations of the real density from these model assumptions are reflected in the residual input signals of RLSC. Accordingly, they are consistently collocated to the output quantities. This issue is not specific for the present RLSC method, but holds generally for all regional gravity field modeling methods that are using the RCR technique. Yang et al. (2018) demonstrate that the use of available lateral density maps does not necessarily improve the gravity modeling results.
The reduced gravity disturbances δg which enter RLSC are presented in Fig. 5 and described in the following. For the terrestrial observations (Fig. 5a), we calculate from Eq. 13 with the SH degrees n XGM2018 ∈ {2, 719} and n dV_ELL_Earth2014 ∈ {720, 2159}. ERTM2160 is applied from SH degree 2160 to its maximum resolution, which equals approximately 250 m. Correspondingly, we have for the airborne observations (Fig. 5b) with n XGM2018 ∈ {2, 719} and n dV_ELL_Earth2014 ∈ {720, 5400}. The reduction to different SH degrees is valid in our case since the gravity signal above SH degree 5400 is negligible at the position of the airborne observations. This is first due to signal attenuation with altitude, which significantly reduces the high-frequency gravity signal at flight height. Secondly, the low-pass filter removes approximately 90% of the still remaining signal above SH degree 5400 (compare to Fig. 3). In Eqs. 27 and 28 we apply globally unbiased models to a very high degree. Accordingly, we can safely assume that they do not have a significant bias in our study area, which is a requirement for calculating the height anomaly output ζ out of RLSC (Eq. 13). The final height anomalies ζ out are calculated by restoring the effect of the corresponding models (XGM2018, dV_ELL_Earth2014 and ERTM2160) with exactly the same SH degrees as in the terrestrial reduction. The subscript 'out' refers to specified output points on the Earth's surface which are a combination of a regular 1 × 1 grid and the GSVS17 points. The topographic heights of the output area are presented in Fig. 6, whereby the highest mountains reach more than 4000 m. The location of GSVS17, a west-east traverse of approximately 350 km through the mountains of Colorado, is added in blue. Note that the location of the grid is centered within the area of terrestrial observations, but significantly shifted in relation to the airborne observations (Fig. 5).
In the RLSC approach, we include the original full covariance information of XGM2018. The SH coefficients of XGM2018 and XGM2016 are very similar, but the corresponding covariance information of XGM2018 has been improved and is even more realistic now. The covariance matrices Cˆlˆl, Cˆsˆl and Cˆsˆs of the reduction model (Eqs. 24,26) are calculated from a model covariance function (MCF) that fits the reduced gravity disturbance input (Fig. 5). The approach is explained and justified in detail by Willberg et al. (2019).
We present our result in agreement with the reference potential value W 0 of the IHRS definition (Ihde et al. 2017) in the mean-tide system. Since output points within the JWG are specified on the Earth's surface, we include the theory of Molodensky (Hofmann-Wellenhof and Moritz 2006) and prefer the calculation of height anomalies instead of geoid heights. Comparisons of the absolute potential values and geoid heights, which are also calculated within the JWG, are not discussed here, but presented in Sánchez et al. ((2020), this issue) and Wang et al. ((2020), this issue). However, we also include the calculation of gravity disturbances at the input positions (with a corresponding restore step) in order to compare them to the original measurements in Sect. 5. Fig. 7 Comparison of RLSC result against available models. a the height anomaly ζ out . b Difference to EGM2008 [ζ out -ζ out (ERTM2160)ζ out (EGM2008)]

Analysis of the RLSC grid
In this and the following section, we present and quantify our height anomaly result ζ out (Eq. 29) related to 'the 1 cm geoid experiment'. We refer to Wang et al. ((2020), this issue) for general conclusions and further results from the JWG. The high-resolution measurements along GSVS17 have not been published yet, therefore we include internal and external validation to verify our results. At first, in Fig. 7, we compare our result against available models to indicate main effects in the result. Fig. 7a shows the height anomaly ζ out from RLSC (without restore step), which gives the improvement from the gravity measurements with respect to the prior models. As expected, we see a high correlation between the output residuals ζ out and the input residuals in Fig. 5. The dominating effects are signals with a spatial extent of approximately 0.1 • to 0.3 • which corresponds to SH degrees between 720 and 2160. In these frequencies, the reduction model consists only of topography-derived gravity information. Thus, we assume that the signals in Fig. 7a are significant improvements due to the input gravity measurements. However, there is also one maximum on the left-hand side with a much larger spatial extent. The terrestrial input residuals (Fig. 5a) especially show a positive bias in this area which seems to be responsible for this maximum. The spatial extent of the red area in Fig. 7a indicates that XGM2018 and the terrestrial gravity observations are not consistent in this area. Figure 7b presents the difference between the height anomaly result ζ out and the corresponding height anomaly derived from EGM2008 (Pavlis et al. 2012). However, in this comparison we do not restoreζ out (ERTM2160) in Eq. 29, because the resolution of ERTM2160 is not included in EGM2008. The presented difference is dominated by features with a spatial extent in the order of 100 km and can be attributed to the difference between EGM2008 and our three data sources: XGM2018, the terrestrial, and airborne gravity observations. Note that the dominating effects from Fig. 7a are not visible in Fig. 7b. We therefore conclude that EGM2008 represents the gravity field between SH degree 720 and 2160 significantly better than the mere topographic gravitational model dV_ELL_Earth2014.
In a second step, we calculate gravity disturbances for the input points of RLSC and compare them to the original observations. We refer to the appropriate difference as output residuals, whereby their presentation in Fig. 8 is limited to the area of the output grid. The output residuals indicate the improvement of the combined RLSC solution with respect to the original gravity measurements. For the terrestrial output residuals in Fig. 8a, we see significant differences in the target area: In areas with topographic heights below 2500 m, e.g., east of − 105% longitude or in the south-west corner, the output residuals are close to 0 mGal. In the mountainous areas of Colorado they often exceed ± 5 mGal. It is common for LSC methods to smooth the input gravity observations in regions with strongly varying gravity signals (e.g., mountains). In our case, this is the result of a homogeneous and isotropic covariance function above the maximum resolution of XGM2018, which consequently leads to high output residuals. However, the very inhomogeneous parts in Fig. 8a indicate that at least some of the output residuals are outliers or measurement errors in the gravity database.
At this point, we do not detect or exclude outliers in RLSC, as there is no beneficial metadata. Furthermore, individual data point selections will inevitably complicate the comparisons with other groups from 'the 1 cm geoid experiment'. However, we acknowledge that an additional, individual data point inspection and corresponding outlier detection could be able to improve the final results. Once more, the area with an offset to XGM2018 stands out in the center of the mid left side, where a large area is dominated by positive output residuals and warm colors. Apparently, RLSC combines the differing data sets XGM2018 and terrestrial observations, thus resulting in residuals to XGM2018 (visible in Fig. 7) and to the terrestrial gravity observations, respectively (visible in Fig. 8). We see a possible reason for the bias in the inhomogeneous distribution of the terrestrial observations which are more often located in mountain valleys. The standard deviation of all terrestrial output residuals is 2.3 mGal, Fig. 8 Output residuals in terms of gravity disturbance, namely original observations−RLSC result for a the terrestrial gravity observations and b the airborne gravity observations and its mean value is 0.3 mGal. The output residuals for the airborne gravity observations, in terms of the original observations, have a standard deviation of 1.6 mGal with a mean value of − 0.1 mGal (Fig. 8b). The dominating effects are again short-wavelength effects in the mountainous regions of Colorado, once again affected from the smoothing effect of LSC. However, this time the output residuals are affected from the included low-pass filter as well. Furthermore, differences between terrestrial and airborne gravity observations are assumed to be a main contributor for the output residuals in Fig. 8. As there are no long wavelength signals visible in Fig. 8b, we conclude that XGM2018 and airborne gravity observations are consistent in the target area.
The residual's standard deviations confirm that the scale of our a priori observation accuracy, σ ter = 3 mGal and σ air = 1 mGal respectively, are generally reasonable. However, the data sets seem to show a spectral dependence, as indicated by much higher residuals in the Colorado mountains. As a result of the equally weighted observation accuracies for a single data type, RLSC produces output accuracies depending on the data distribution and a priori accuracies only. Accordingly, the estimated accuracies cannot adequately represent the problems of suspicious areas that we see from the residuals. Figure 9 shows the resulting 3σ confidence level for height anomalies in the target area which is derived from the error covariance matrix E ss in Eq. 26. The values depend on the accuracy assumptions for the gravity observations, their point distribution and the covariance information of XGM2018. The 3σ confidence level varies from approximately 1 cm in areas with very dense terrestrial observations to more than 5 cm in the very north where the solution is not supported by airborne observations (cf. Fig. 5). The availability of accuracy estimates is a main advantage of LSC approaches. However, in this case an iterative RLSC approach would be necessary to derive more realistic output accuracies. For example, we could use the residuals in Fig. 8 to derive individual input In summary, the terrestrial gravity observations show notable differences to XGM2018 and airborne gravity observations in the highly mountainous regions of Colorado, but a generally good consistency everywhere else. We assume that the height anomaly result ζ out could be further improved by an outlier detection, or far-reaching data inspection in general. Additionally, an iterative RLSC approach could derive accuracy assumptions from the output residuals. This would result in a more selective weighting between the data sets and an improved accuracy estimation in general. However, we see these aspects as beyond the scope of this paper. Our focus is, first, on the derivation of a consistent methodology for lowpass filtering within RLSC, and second on the generation of a comparable result within 'the 1 cm geoid experiment'.

Comparison along GSVS17
The height anomalies along GSVS17 have been computed by 13 different groups and are presented with different colors in Fig. 10a. The values at the 223 benchmarks are given as difference to the joint mean value of each individual point. As long as the real reference values along GSVS17 are not yet submissions with an individually corrected mean value, whereby our estimated confidence interval is additionally shown in gray color available, we assume that the mean value of 13 independent calculation methods is significantly better than an individual solution. This notion is, for example, used analogously in meteorological literature (Evans et al. 2000;Ebert 2001). Therefore, we interpret the common mean value as a reference which is represented by zero in Fig. 10. However, it should be noted that this mean value is not necessarily without systematic effects, as they could be introduced by the input gravity observations, and thus reflected in all solutions. During several phases of the JWG the offset between different solutions was reduced , this issue), mainly by the adaptation of common standards, e.g., zerodegree term and tide system. The remaining difference in the mean value might result from differences in the processing strategies, the topographic reduction or the individual data handling, and cannot be solved by this paper. At the present stage, the absolute mean offsets of the different curves in Fig. 10a range from 0.1 to 2.7 cm, and the standard deviations from 7.5 to 2.4 cm. Our solution, which is highlighted by a black color, has a mean value of − 1.0 cm, and with 7.5 mm the lowest standard deviation among all results. As a result, the variations in our solution are significantly reduced in comparison to most others (Fig. 10a).
In Fig. 10b, the remaining issue of the mean offset between different height anomaly solutions is ignored and individually corrected. Furthermore, we selected only seven out of the 13 solutions which are most consistent to the previously introduced mean value. Our height anomaly result (black line) is the only solution that varies only within ±2 cm in relation to the joint mean value. The estimated 3σ confidence interval of our solution is additionally overlaid with gray color. We conclude that several of the other results are mainly covered within this confidence level, in spite of the fact that our provided standard deviation is probably too optimistic in the mountainous parts of Colorado (see Sect. 5). Even when we reduce a mean offset and consider only the solutions in Fig. 10b, there are remaining variations of some centimeters between the different height anomaly results. Note that similar approaches might benefit from the comparison to a mean value as they are more likely to end up close together. However, our approach is the only RLSC approach, and the only one-step LSC method among the contributing groups.
One of the main reason for the remaining variations in Fig. 10 are the differences between terrestrial and airborne observations in the area. In order to quantify the spread resulting from different weighting between terrestrial and airborne observations, a series of additional test computations have been performed. In Fig. 11, the black line is again our height anomaly result ζ out in reference to the mean, while other curves represent alternative solutions. We exemplary include the following four cases: (1) using only terrestrial observations and disregarding airborne data in blue, (2) using only airborne observations in red, (3) using a different airborne input accuracy in green, and (4) a solution without low-pass filter of airborne data in purple. In general, significant differences show up, especially between the terrestrial-only and the airborne-only result. In an extreme case (around point ID 110) two airborne flight lines are positioned above two opposing mountain flanks, while GSVS17 runs through the in-between valley with a distance of approximately 5 km to both of the flight lines. As a result, the two corresponding solutions differ by more than 8 cm, obviously measuring a different gravity signal. We conclude from the height anomaly results shown in Fig. 11 that rather large differ-ences of several centimeters can be caused by the weighting of the individual data sets.
As an example, the green curve in Fig. 11 is calculated equally to the black one, but uses a different input accuracy for the airborne observations, namely σ air = 1.5 mGal instead of 1 mGal. As a result, RLSC gives a lower relative weight to the airborne observations, which leads to a shift of the combined solution towards the terrestrial-only solution. An opposite effect is visible for the purple curve. It is a combined solution that includes unfiltered airborne observations with an 8 Hz sampling frequency and an accuracy assumption of σ air = 2 mGal since the high-frequency noise is still included. The significantly increased number of airborne observations, and the fact that they are assumed to be uncorrelated, amplifies their relative weight in RLSC. Accordingly, this solution is shifted towards the airborne-only solution (red curve). This could be an indication that one of the main reason for the height anomaly differences in Fig. 10a is the various treatment of terrestrial and airborne input data. Therefore, it is not only related to the different regional gravity modeling approaches used in this inter-comparison exercise.
In summary, our height anomaly ζ out performs very well in a comparison among the JWG results and shows the smallest variations in regard to the joint mean value. We highlight differences in the available data sources, and conclude that their relative weighting will be one of the main drivers of the final performance. However, we emphasize that all data issues mentioned in 'the 1 cm geoid experiment' are not a problem of our gravity modeling approach, but rather mere inconsistencies in the available data sources. In order to solve these issues, a detailed data inspection of the original data sources would have to be done, which would go far beyond the scope of this study. However, solving some of the data issues could become easier once the high-quality observations from the GSVS17 project are released. At the present stage, we conclude that different scenarios of RLSC show consistent results in a simple comparison along the GSVS line. Furthermore, we demonstrate that our method provides the necessary flexibility for an adjustment in case a data set provides either benefits or problems in a target area.

Conclusion and outlook
In this paper, we use a low-pass filter for reducing highfrequency noise in airborne gravity observations from the GRAV-D project. Accordingly, we derive a novel concept in order to include the resulting correlations to regional gravity field modeling. In RLSC, the functional and stochastic model is adapted, whereby filtered observations and filtered covariance matrices are treated in a consistent manner. The approach is verified with a combination of filtered (airborne) Fig. 11 Comparison of different internal height anomaly results along GSVS17. Solutions vary according to their relative weighting between terrestrial and airborne observation, and show differences of up to 8 cm and unfiltered (terrestrial) observations in the frame of 'the 1 cm geoid experiment'.
In the target area of Colorado, our calculated height anomaly grid displays significant improvements compared to already existing gravity models. They are visible as long wavelength differences to the EGM2008, and short to medium-wavelength deviations from the a priori topographic gravity model. The calculated output residuals for the gravity disturbance indicate a very good consistency of terrestrial and airborne data in areas with topographic heights below 2500 m, but reveal some issues in the highly mountainous regions of the Rocky Mountains. We even identified one region where the terrestrial observations differ significantly from the long wavelength part of XGM2018.
Within 'the 1 cm geoid experiment' the RLSC method performs very well in a comparison among 13 independent height anomaly results along the GSVS17 benchmarks. Since the actual results of GSVS17 are not yet available, the reference value which is used to evaluate the performance of our solution, is calculated from a common mean value. In general, the final height anomaly results deviate from the reference by a few centimeters. With a standard deviation of only 7.5 mm our solution shows the smallest variations with respect to the reference. Furthermore, it is the only result that always stays within ± 2 cm (in case the individual offsets are reduced).
Moreover, we analyzed the impact of the systematic differences between the available data sources on the final result. In a comparison along GSVS17, we tested four further data and processing scenarios, and could show that the relative weighting among the input data types can cause differences in the order of several centimeters in the height anomaly results. Consequently, RLSC can be individually adapted in the case that there is additional information about the quality of the gravity observations. More generally, we conclude that a significant part of the differences in 'the 1 cm geoid experiment' might be related to a different treatment and relative weighting of the input data, and not only to various regional gravity field modeling methods. In this regard, it might be useful to check the internal weighting between our data sets again, as soon as the high-quality measurements along GSVS17 are available to the public.
One of the main advantages of the statistical method of (R)LSC is that it also provides variance-covariance information of the resulting quantities. The availability of realistic error estimates will be crucial for the realization of the IHRS. Due to restricted availability and heterogeneous quality of ground gravity data, it is to be expected that a similar accuracy of absolute potential values at the IHRS stations cannot be achieved. In this case, the provision of realistic error estimates in addition to the potential values themselves will be very important. The estimated standard deviations of height anomalies are, in our current solution, mainly dominated by the distribution of the gravity observations as well as the relative weighting of the input data sets. However, a constant weight of the two data sets, terrestrial and airborne observations, has been included. In an iterative RLSC approach, the post-fit residuals could be used to modify the a priori accuracy of the input data, and thus the relative weighting scheme. By this, the identified systematic differences between terrestrial and airborne data in certain areas could also be taken into account. Consequently, this would result in an even more realistic stochastic modeling and further improved error estimates for the output quantities. Since the main goal was to inter-compare the results with other study teams, this has not been applied in the present solution but will be part of further studies.
Summarizing, 'the 1 cm geoid experiment' helps for the scheduled IHRS definition and realization because it provides a meaningful accuracy benchmark in the case of a good data distribution in very difficult terrain. Within this JWG all contributing groups have proven their capability to calculate the height anomaly with an accuracy of some centimeters.
In the end, IHRS should be defined by absolute potential values, but their derivation from the disturbing potential T , which is also needed for the height anomaly, is straightforward. Sánchez et al. ((2020), this issue) propose a globally distributed computation of potential values at IHRS reference stations due to common data restrictions. Thereby, for example, all institutions with a specific performance in terms of a reference could be allowed to contribute to the calculation of IHRS stations with their own method. Consequently, it will be extremely important to analyze, how much of the remaining deviations are caused by differences in the regional gravity modeling methods, and how much by the different treatment and relative weighting of input data.
The benefit of including covariance information from high-resolution GGMs in RLSC and the resulting advantage in terms of the calculation of IHRS reference stations is already explained in Willberg et al. (2019). With the consistent handling of airborne gravity data and by including a covariance propagation of the filter behavior, we can further enhance the stochastic modeling, which, as discussed above, will be very beneficial for the IHRS realization.