The RTM harmonic correction revisited

In this paper, we derive improved expressions for the harmonic correction to gravity and, for the first time, expressions for the harmonic correction to potential and height anomaly. They need to be applied at stations buried inside the masses to transform internal values into harmonically downward continued values, which are then input to local quasi-geoid modelling using least-squares collocation or least-squares techniques in combination with the remove-compute-restore approach. Harmonic corrections to potential and height anomaly were assumed to be negligible so far resulting in yet unknown quasi-geoid model errors. The improved expressions for the harmonic correction to gravity, and the new expressions for the harmonic correction to potential and height anomaly are used to quantify the approximation errors of the commonly used harmonic correction to gravity and to quantify the magnitude of the harmonic correction to potential and height anomaly. This is done for two test areas with different topographic regimes. One comprises parts of Norway and the North Atlantic where the presence of deep, long, and narrow fjords suggest extreme values for the harmonic correction to potential and height anomaly and corresponding large errors of the commonly used approximation of the harmonic correction to gravity. The other one is located in the Auvergne test area with a moderate topography comprising both flat and hilly areas and therefore may be representative for many areas around the world. For both test areas, two RTM surfaces with different smoothness are computed simulating the use of a medium-resolution and an ultra-high-resolution reference gravity field, respectively. We show that the errors of the commonly used harmonic correction to gravity may be as large as the harmonic correction itself and attain peak values in areas of strong topographic variations of about 100 mGal. Moreover, we show that this correction may introduce long-wavelength biases in the computed quasi-geoid model. Furthermore, we show that the harmonic correction to height anomaly can attain values on the order of a decimetre at some points. Overall, however, the harmonic correction to height anomaly needs to be applied only in areas of strong topographic variations. In flat or hilly areas, it is mostly smaller than one centimetre. Finally, we show that the harmonic corrections increase with increasing smoothness of the RTM surface, which suggests to use a RTM surface with a spatial resolution comparable to the finest scales which can be resolved by the data rather than depending on the resolution of the global geopotential model used to reduce the data.

cable, the bathymetry) appeared to be necessary to obtain an accurate local model. The former ensures that the long wavelengths have less energy, which makes it possible to ignore data far outside the area of interest when computing the local quasi-geoid model. The latter ensures that the energy of short wavelengths (ideally those which cannot be modelled) is significantly reduced, and the residual field, being harmonic in a larger domain, is smoother than the original one.
There are several types of topographic reduction, depending on the method used in local quasi-geoid modelling, see Tziavos and Sideris (2013) for a review. The subject of this study is the so-called residual terrain model (RTM) reduction, which has been introduced in (Forsberg and Tscherning 1981). This type of topographic reduction is part of the remove-compute-restore approach and frequently used in the framework of Molodensky's theory of the quasi-geoid (Denker 2013;Tziavos and Sideris 2013). Given the fact that in local quasi-geoid modelling, the data are always reduced for the contribution of a GGM, the aim of the RTM correction is to reduce the data for the gravitational signal of the masses, which are not captured by the GGM. To do so, the real Earth's topography (including bathymetry if applicable) is replaced by a smoother topography (the RTM topography), and the data are reduced for the gravitational effect of the mass deficit between the two topographies. In the past, this idea was implemented in different ways (cf., Forsberg and Tscherning 1981;Forsberg 1984;Denker and Wenzel 1986;Vermeer and Forsberg 1992;Hirt et al. 2014); see also the review paper (Tziavos and Sideris 2013).
In this study, we follow the implementation in (Forsberg and Tscherning 1981). The body with RTM topography is referred to as the RTM-reduced earth and its surface is referred to as the RTM surface (cf. Fig. 1). A wellknown conceptual problem of this implementation arises for data points located below the RTM surface, a situation frequently encountered in or close to areas of strong topographic variations. After RTM reduction, the gravity field functionals neither represent boundary values nor values inside the harmonic domain of the gravitational potential of the RTMreduced earth. To address this conceptual problem, Forsberg and Tscherning (1981) suggested an additional correction to be applied to gravity anomalies at those stations, which they refer to as the harmonic correction. After RTM and harmonic correction have been applied, the reduced gravity anomalies are consistent with a harmonically downward continued gravitational potential of the RTM-reduced earth. The downward continuation is not done strictly, but comprises several approximations. The most important one is that for every buried data point the action of the masses encompassed in the RTM body is modelled as a Bouguer plate of a thickness equal to the vertical distance of the data point to the RTM surface. That is, to every data point a planar Bouguer plate is assigned with its own vertical position and thickness. The advantage of this simplification is that the harmonic correction can be computed easily (cf. Sect. 2). The downside is that errors in the reduced gravity anomalies are introduced, which among others, depend on the deviation of the RTM surface from a planar surface, in particular in some neighbourhood of the data point. The deviation depends on how the RTM surface has been constructed from a high-resolution model of the topography and bathymetry. A common approach is to apply a low-pass filter, and to relate the cut-off frequency of this filter to the maximum spherical harmonic degree of the GGM. Then, the higher the maximum degree is, the more the RTM surface differs from the surface of a Bouguer plate. Errors in the reduced gravity anomalies are also introduced if the Bouguer plate is replaced by a finite Bouguer plate, a spherical shell or a spherical cap. (e.g., Kadlec 2011). Forsberg and Tscherning (1981) do not apply a harmonic correction to disturbing potential or height anomaly (see also Forsberg 1984). They argue that if the RTM surface is smooth and does not differ much from a Bouguer plate in the neighbourhood of the station, the harmonic correction for potential and height anomaly are close to zero, and may be neglected.
The main objective of this paper is twofold. First, for stations located below the RTM surface, we derive new formulas for harmonically downward continued gravity and potential values referring to the RTM-reduced earth bounded by the RTM surface. Second, we use the new formulas and investigate (i) the approximation error of the harmonic correction to gravity as suggested in (Forsberg and Tscherning 1981), and (ii) the magnitude of the harmonic correction to potential and height anomaly, which were neglected until now. We do this numerically for two areas with different topographic regimes.
The paper is organised as follows. In Sect. 2, we briefly summarise and elaborate on the harmonic correction to gravity in (Forsberg and Tscherning 1981). In Sect. 3, we derive the new formulas and extract from them new expressions for the harmonic correction to gravity and potential and compare them with the harmonic correction in (Forsberg and Tscherning 1981). In Sect. 4, we describe the results of experiments, which aim at a verification of the formulas derived in Sect. 3, and a quantification of (i) the approximation error of the harmonic correction to gravity in (Forsberg and Tscherning 1981) and (ii) the magnitude of the harmonic corrections to potential and height anomaly. In Sect. 5, we provide a summary and some concluding remarks.

Harmonic correction of Forsberg and Tscherning
The RTM correction in (Forsberg and Tscherning 1981) reduces an Earth gravity field functional to the corresponding functional of a RTM-reduced earth, which is bounded by the RTM surface, s. A point P ∈ S located above the RTM surface hangs in free air whereas a point P ∈ S located below the RTM surface is buried inside the masses of the RTM-reduced earth (cf. Fig. 1). Hence, a reduced gravity field functional at a point P ∈ S located below the RTM surface represents a functional inside the masses of the RTM-reduced earth. The harmonic correction of (Forsberg and Tscherning 1981) aims at a transformation of this internal gravity field functional into a harmonically downward continued gravity field functional. Forsberg and Tscherning (1981) apply this transformation only to gravity. Let P ∈ S be a point on the Earth's surface (i.e., on the topography or the sea surface); let Q be the intersection of the ellipsoidal normal through P with the RTM surface, i.e., Q ∈ s. As we assume that P is located below the RTM surface, the ellipsoidal height difference h = h Q − h P is positive. The transformation from an internal gravity value to a harmonically downward continued gravity value according to Forsberg and Tscherning (1981) is done in four steps (e.g., Gerlach 2003, p. 71).
Step 1 Add 2π Gρ 0 h to the gravity at P, g(P). This corresponds to the removal of a Bouguer plate of thickness h and mass density ρ 0 , which is located on top of point P.
Step 2 Use a model of the free-air gravity gradient along the ellipsoidal normal through P, and upward continue gravity from P ∈ S to Q ∈ s. This provides gravity at Q, g(Q).
Step 3 Restore the Bouguer plate of step 1; this increases g(Q) of step 2 with 2π Gρ 0 h. Note that steps 1 to 3 are comparable to the Poincaré and Prey reduction used in Helmert's theory ( Heiskanen and Moritz 1967, Sect. 4.3).
Step 4 Downward-continue g(Q) of step 3 using the free-air gravity gradient of step 2. This provides a gravity value at P, which is consistent with a harmonically downward continued potential of the RTM-reduced earth with boundary s.
After step 4, gravity g(P) has increased with Equation (1) is the harmonic correction to gravity in (Forsberg and Tscherning 1981). Note that although it is called a correction, it is added to the RTM-reduced gravity value to obtain a harmonically downward continued gravity value, whereas the RTM correction is subtracted from measured gravity. In this paper, we follow the sign convention in (Forsberg and Tscherning 1981) for the harmonic correction. The increase of Eq. (1) is twice the effect of a Bouguer plate of thickness h and density ρ 0 . The same result is obtained if in step 3 the Bouguer plate is not restored but condensed on an infinite plane through point P. This interpretation is used in (Forsberg and Tscherning 1981). Applying the same procedure to a height anomaly is meaningless as the gravitational potential of a Bouguer plate is not defined. However, Forsberg and Tscherning (1981) argue that if the RTM surface is smooth enough in some neighbourhood of the data point, the harmonic correction to potential (likewise to height anomaly) at this point is very small and may be neglected. Ignoring the harmonic correction to potential implies that no difference is made between the potential inside the masses and a harmonically downward continued potential. Note that Eq.
(1) still applies if the Bouguer plate is replaced by a spherical shell. This follows directly from the formulas derived in Kadlec (2011) and is a result of potential theory (Martinec 1998). However, in case of a spherical shell, the harmonic correction to height anomaly is well-defined and given by where γ is normal gravity at P (cf. Kadlec 2011). Assuming ρ 0 = 2670 kg/m 3 and h = 500 m, this is about 2.8 cm, i.e., not negligible in cm-accuracy quasi-geoid modelling. The effect is below 1 cm if | h| < 300 m. This simple example already indicates that depending on the situation, the harmonic correction to height anomaly may not be neglected.
The harmonic correction to gravity in (Forsberg and Tscherning 1981) is an approximation, and to our knowledge, the approximation error has not been investigated yet. A major conceptual weakness of this approximation is that the vertical position of the Bouguer plate differs per data point. Even if two data points are just a few kilometres apart, the vertical positions of the two Bouguer plates may differ by tens or hundreds of metres. This introduces inconsistencies in the reduced gravity values, in particular between neighboured data points.
There are several attempts in literature to address the conceptual weakness of the harmonic correction in (Forsberg and Tscherning 1981). For instance, Kadlec (2011) suggested to replace the (infinite) Bouguer plate by a finite one or a finite spherical shell (i.e., a spherical cap), and writes the RTM correction as the difference between the complete effect of the topography and the complete effect of the RTM topography. However, this does not provide gravity field functionals consistent with a harmonically downward continued potential field. Therefore, this idea should not be used in the context of local quasi-geoid modelling. Omang et al. (2012) suggested to upward continue inside the masses observed gravity field functionals to the RTM surface ignoring the horizontal components of the gravity gradient vector. After RTM reduction, the gravitational potential of the RTM-reduced earth fulfils Laplace's equation, i.e., V (P) = −4π Gρ 0 , and assuming vanishing horizontal gravity gradients, it is V zz = −4π Gρ 0 . Therefore, the upward continuation to the RTM surface changes gravity by −4π Gρ 0 h. In this way, they obtain gravity field functionals at points on the RTM surface, and avoid the problem of a harmonic correction. The main disadvantage of this approach is that errors in the vertical gradient of gravity, specially due to variations of density, directly propagate into the upward continued gravity anomalies proportional to the height difference h.

The complete RTM correction
Here, we follow a different approach for data points located below the RTM surface. The approach is applicable to all methods used in local quasi-geoid modelling which do not require data located on the boundary of the harmonic domain. The approach does not move terrestrial data to the RTM surface. This is important, because moving data from the Earth's surface (land topography or sea surface) to the RTM surface may introduce significant errors in areas where the distance to the RTM surface is large (e.g., the Norwegian fjords) due to a lack of knowledge of the external potential field. This applies in particular to gravity anomalies due to uncertainties in the free-air vertical gravity gradient. Furthermore, the approach provides harmonically downward continued data referring to a RTM-reduced earth which is bounded by the RTM surface like the harmonic correction in (Forsberg and Tscherning 1981) does, but with a much higher accuracy and not limited to gravity anomalies. This aspect is very important if gravity data are combined with other gravity field functionals in local quasi-geoid modelling as only then the reduced datasets are consistent in the sense that they all refer to a harmonically downward continued field.
Intrinsically, the new equations not only comprise the RTM corrections but also improved harmonic corrections. We refer to them as the complete RTM corrections. The equations are derived for two types of gravity field functionals, i.e., for gravity g and potential W . Dividing the equation for W by the normal gravity at the telluroid point provides the equation for height anomaly. The derivation of the corresponding equations for other gravity field functionals is out of the scope of this study, but straightforward. We find it appropriate to note that for all other points (i.e., points located above the RTM surface), the complete RTM correction is identical to the RTM correction in (Forsberg and Tscherning 1981).
We assume that the data point P is located on the land topography or the sea surface. As before, a point Q is the intersection of the ellipsoidal normal through P with the RTM surface. denotes the volume between the Earth's surface (i.e., the topography on land or the bathymetry at sea) and the RTM surface. The part of the volume that is located above the RTM surface is denoted + ; likewise, the volume below the RTM surface is denoted − . The RTM-reduced earth is mass-free in volume + and is filled with masses of constant density ρ 0 in volume − , where ρ 0 is commonly set equal to a mean crustal density (cf. Figs. 1, 2). δV + RTM and δg + RTM denote the gravitational potential and gravitational strength, respectively, of the masses inside the volume + . Likewise, δV − RTM and δg − RTM denote the gravitational potential and gravitational strength, respectively, of the masses added to the volume − .
Note that when computing δg + RTM , δg − RTM , δV + RTM , and δV − RTM we need to distinguish between the five cases depicted in Fig. 2 and listed in Table 1. Three of them require to compute two tesseroids per grid cell (Heck and Seitz 2007;Grombein et al. 2013).

Complete RTM correction for points located below the RTM surface
Let δg corr (P) and δW corr (P) be the complete RTM correction to measured gravity and potential, respectively. Then, we define the reduced gravity g red (P) and reduced potential W red (P) as Note that g red (P) and W red (P) refer to a RTM-reduced earth, which is mass-free outside the RTM surface; moreover, they represent harmonically downward continued potential field functionals at points P ∈ S located below the RTM surface. Whenever convenient, the value of a function f at a point P is either written as f (P) or as f P . All equations are provided with error terms. The factors appearing in the error terms are derived for an isotropic potential field for which applies The complete RTM correction requires four steps to be taken: Step 1 Move the masses inside + to infinity and correct gravity and potential for the gravitational effect this has. If the effect on gravity is denoted δg + RTM (P) and the effect on potential is denoted δV + RTM (P), we obtain reduced values of gravity, g + (P), and potential, W + (P), respectively, i.e., Note that all evaluation points P are located on the Earth's surface (land topography or sea surface), i.e., P ∈ S (cf. Figs. 1, 2).
Step 2 For points P ∈ S located below the RTM surface (i.e., points for which h = h Q − h P > 0, where h is the ellipsoidal height), upward continue gravity and potential obtained in step 1 to the corresponding points Q ∈ s; data at all other points P ∈ S are left unchanged: where R is the mean radius of the Earth and ε = h R . Terms of the order of O(ε 2 ) and O(ε 3 ) are neglected in the expansion of gravity and potential, respectively. Note that ∂ g + ∂h in Eq. (8) and ∂ 2 W + ∂h 2 in Eq. (9) have step discontinuities across the surface bounding the masses. Therefore, it is appropriate , ∂ W + ∂h P and ∂ 2 W + ∂h 2 P are freeair gradients which apply to a modified Earth which has no masses in + , i.e., where P is a point on the ellipsoidal normal through P with ellipsoidal height h P + δ. Using the difference quotient, Eqs. (10) and (11) may be approximated by respectively.
Step 3 Fill the volume − with mass so that the mass density equals ρ 0 , and correct gravity and potential for the gravitational effect this has (cf. Figs. 1, 2). For data that have been upward continued in step 2, the gravitational effect has to be computed at the points Q ∈ s. If δg − RTM (Q) is the effect on gravity and δV − RTM (Q) the effect on potential, we obtain reduced values of gravity and potential, respectively, i.e., Note that for all other points, the reduced values of gravity and potential refer to points P ∈ S.
Step 4 Perform a harmonic downward continuation of the data referring to points Q ∈ s (i.e., data that have been upward continued in step 2) to their original locations P ∈ S on the Earth's surface (i.e., land topography or sea surface) along the ellipsoidal normal through P: The gradients in the last two equation are free-air gradients which apply to a RTM-reduced earth bounded by the RTM surface, i.e., and Q is a point on the ellipsoidal normal through Q with ellipsoidal height h Q + δ. Hence, for points P located below the RTM surface, g red (P) and W red (P) represent harmonically downward continued quantities.   (14) into Eq. (16) yields Likewise, inserting Eqs. (7), (9), and (15) into Eq. (17) yields In Eq. (21), the term proportional to h can be simplified.
To do so, we take the partial derivative with respect to h of Eqs. (14) and (8), respectively, and obtain, Inserting the last equation into Eq. (21) gives From Eq. (14), which allows to simplify the error term of Eq. (24) and provides the final form of Eq. (21): The term in brackets on the right-hand side of Eq. (26) is the complete RTM correction to gravity as defined in Eq. (3). Ignoring the error term provides the approximation Equation (22) can be simplified. First, taking the first and second partial derivative of Eq. (15) with respect to h, we find Inserting this result into Eq. (22) gives A Taylor series expansion of W + about P and Q, respectively, shows that in Eq. (29) the last term before the error term is of the order of 4W + (Q)O(ε 4 ), i.e., it can be neglected. Furthermore, the error term in Eq. (29) can be simplified. From Eq. (15), it follows that hence, which provides the final from the Eq. (22): The term in brackets on the right-hand side of Eq. (32) is the complete RTM correction to potential as defined in Eq. (4).
Ignoring the error term provides the approximation We may achieve the same results within error in just three steps. This is shown in "Appendix A".

Complete RTM correction for arbitrary data points
So far, we derived new formulas for the complete RTM correction for data points located below the RTM surface. For all other data points, the well-known RTM correction in (Forsberg and Tscherning 1981) applies. Hence, we obtain the following formulas for arbitrary data points P ∈ S: , P is located above the RTM surface and therefore in free air after mass redistribution), the reduced value of gravity at this point is The term in brackets is identical to the RTM correction given in Forsberg and Tscherning (1981) (see also Forsberg and Tscherning (1997)). At a point P ∈ S with h = h Q − h P > 0 (i.e., P is located inside the masses of the RTM-reduced earth), the reduced value of gravity Norway test area: harmonic correction to gravity, δg harm (P) (Eq. (37)), the difference between δg harm (P) and the harmonic correction in (Forsberg and Tscherning 1981), δg harm (P) − δg FT harm (P), and the harmonic correction to height anomaly, δζ harm (P) (based on Eq. (41)). The areas where topography or sea surface is located above the RTM surface are shown in white We re-write Eq. (35) as The first two terms on the right-hand side of Eq. (36) represent a gravity value inside the masses of the RTMreduced earth. Hence, the last term in brackets reduces this internal gravity value into a harmonically downward continued gravity value, i.e., it is the new formula for the harmonic correction to gravity: Equation (37) the reduced value of the potential is If P ∈ S with h = h Q −h P > 0, the improved expression for the reduced potential up to terms δV − Formally, Eq. (39) may be re-written as where the second term in brackets is the new formula for the harmonic correction to potential, i.e., Equation (41) (37) and (41) The former is much larger than the accuracy of today's gravity anomaly datasets; the latter is below the accuracy level achievable today in local gravity field modelling.

Approximation error of the harmonic correction of Forsberg and Tscherning
The complete RTM correction according to Forsberg and Tscherning (1981), δg FT corr (P) and δW FT corr (P), respectively, is the sum of RTM correction and harmonic correction, with the understanding that the latter is only applied to gravity at points located below the RTM surface and is assumed to be zero for potential, i.e., where h = h Q − h P , and δg FT harm (P) is the harmonic correction of Eq. (1) proposed in (Forsberg and Tscherning 1981).
If h > 0, there is a difference between the complete RTM correction δg corr (P) of Eq. (27) and δg FT corr (P) of Eq. (42). Assuming that the approximation error of δg corr (P) is much smaller than that of δg FT corr (P) (this is confirmed by the numerical results of Sect. 4), we may refer to the difference as the approximation error of the harmonic correction to gravity in (Forsberg and Tscherning 1981): Likewise, if h > 0, there is a difference between the complete RTM correction to potential, δW corr (P) of Eq. (33) and δW FT corr (P) of Eq. (43). This difference is equal to the new formula for the harmonic correction to potential as (Forsberg and Tscherning 1981) assumed that the harmonic correction to potential is negligible for smooth RTM surfaces:

Experiments
The experiments are divided into two parts. First, semianalytical investigations are carried out to clarify the nature of the harmonic correction and to verify the equations derived in Sect. 3. Second, numerical investigations are carried out in two regions with the goal to quantify the approximation error of the harmonic correction to gravity and potential in Forsberg and Tscherning (1981). Here, we present the results of the numerical experiments; the semi-analytical experiments are presented in "Appendix B". Numerical experiments were done for two regions with different topographic regimes. The first region is located between 1 • E-10 • E and 56 • N-63 • N and comprises parts of Norway and the Atlantic ocean. Here, we expect significant approximation errors of the harmonic corrections to gravity and potential in Forsberg and Tscherning (1981) as height differences between the Earth's surface and the RTM surface are large and change quickly, in particular around the fjords. The second region is located in the centre of France between 1 • W-7 • E and 43 • N-49 • N. This area is referred to as the Auvergne test area as it has been used in the past to compare different methods of local quasi-geoid modelling (cf. Duquenne 2006). It is a moderate mountainous area with heights varying from less than 150 m in the north-west where the Bassin de Paris starts to 1886 m at the Puy de Sancy in the middle south. In this test area, we expect smaller approximation errors of the harmonic correction to gravity in (Forsberg and Tscherning 1981) and smaller harmonic corrections to potential and height anomaly.
Each RTM surface is computed by applying a low-pass filter to the high-resolution digital topographic and bathymetric model (DTM). We use a spherical Gaussian low-pass filter, which is truncated at a spherical distance ψ 0 . The spherical distance ψ 0 is also the distance at which the Gaussian filter has dropped to half of it maximum value. We refer to ψ 0 as the radius of the Gaussian filter. We relate the choice of ψ 0 to the maximum degree N of the GGM which is used to reduce the data for long-wavelength signals, as ψ 0 = 180 • N , i.e., ψ 0 is set equal to the half-wavelength resolution of the GGM. This filter is in fact a weighted moving average filter, where the size of the window is equal to the half-wavelength resolution of the GGM and the weights are determined by the Gaussian. The filter operation is implemented in the frequency domain using a 2D planar FFT.
For each test area, two RTM surfaces with different smoothness were computed. The first one assumes that the data were reduced for the contribution of an ultra-highresolution GGM complete to degree 2160 (e.g., EGM2008, Pavlis et al. 2012). The corresponding radius of the spherical Gaussian filter is ψ 0 = 5 . The second one assumes that the GGM is complete to degree 300 (e.g., Pail et al. 2011;Bruinsma et al. 2014;Förste et al. 2019); the corresponding radius of the spherical Gaussian filter is ψ 0 = 36 . We refer to the two RTM surfaces as RTM 5 and RTM 36 , respectively.
Applying a low-pass filter to a local DTM causes edge effects along the boundaries of the DTM due to missing data. Therefore, RTM effects were always computed over a spherical patch which was smaller than the area of the DTM by one degree in all four directions.
In Sects. 4.2 and 4.4, we focus on the harmonic correction to gravity, potential, and height anomaly. An analysis of the individual contributors to the complete RTM correction to gravity and potential is provided in "Appendix C" for the Norway test area and in "Appendix D" for the Auvergne test area.

Norway test area
The DTM for the Norway test area is based on EuroDEM (Eurogeographics 2008) and has a half-wavelength resolution of 2 × 2 . A characteristic feature of this area are the long, deep and narrow fjords. Hence, we can expect significant positive and negative height differences between the DTM and each of the two RTM surfaces. These height differences are referred to as RTM tesseroid heights (cf. Fig. 3). The RTM tesseroid heights range from −1404 m to 1036 m for RTM 5 and from −1825 m to 1252 m for RTM 36 (cf. Table 2); the corresponding RMS values are 139 m and 204 m, respectively. In the vicinity of extreme tesseroid heights, we can expect large approximation errors of the harmonic correc-tion to gravity in (Forsberg and Tscherning 1981) and large harmonic corrections to potential and height anomaly.

Norway test area: harmonic correction
In the Norway test area, the harmonic correction was evaluated at the 800,000 nodes of an equal-angular 18 ×18 grid, which formed a subgrid of the DEM. About 32% (RTM 5 ) and 29% (RTM 36 ) of the evaluation points were located below the RTM surface. Figure 4 shows a geographic rendition of the improved harmonic correction to gravity and height anomaly, and the approximation error of the harmonic correction to gravity in (Forsberg and Tscherning 1981). It appears that their distributions are asymmetric with a one-sided long tail, which is the reason why we prefer percentiles and median above RMS, mean, and standard deviation in the summary statistics of Table 3.
Referring to Table 3, there are a few observations worth to be mentioned. First, the harmonic correction to gravity of Eq. (37) can be negative unlike the correction in Forsberg and Tscherning (1981). The latter is always non-negative, as it assumes that all masses of the RTM-reduced earth inside the volume Q − are located above the level of the evaluation point. In reality, however, some masses in − may be located below that level. Overall, this reduces the magnitude of the harmonic correction, but may also lead to negative values at particular evaluation points in extreme situations (e.g., if the RTM surface has lows in the neighbourhood of the evaluation point). In the Norway test area this only happened for the smooth RTM surface, RTM 36 , at just 0.03% of the evaluation points. No negative values were noticed for the rougher RTM 5 surface. We explain this by the fact that then the volume − is smaller so is the amount of mass added to this volume, and the effect of masses located below the level of the evaluation point and close-by is very small and not enough to make the harmonic correction negative.
Second, it is an over-simplification to state that the smoother the RTM surface the smaller the harmonic correction to gravity. For the Norway test area, all statistics of Table 3 point to a larger harmonic correction for the (smooth) RTM 36 surface, compared to the (rough) RTM 5 surface. This can be explained by the fact that for a smooth RTM surface, the volume − increases so does the amount of mass added to that volume resulting in larger harmonic corrections to gravity.
Third, the approximation error of the harmonic correction to gravity in Forsberg and Tscherning (1981) is very large in the Norway test area. It ranges from −126.9 mGal to 13.1 mGal (RTM 5 ) and −284.8 mGal to 191.5 mGal (RTM 36 ). Hence, the range of the approximation error of the harmonic correction to gravity in Forsberg and Tscherning (1981) is comparable to the range of the harmonic correction of Eq. (37). Most of the approximation errors are negative (about 81% for RTM 5 and 77% for RTM 36 ), i.e., the harmonic correction to gravity in Forsberg and Tscherning (1981) is often too large. We explain this by the fact that Forsberg and Tscherning (1981) assumed that all masses of the RTM-reduced earth are located above the level of the evaluation point, whereas in reality there are also masses located below that level. Fourth, the harmonic correction to height anomaly (likewise to potential) is non-negative for both RTM 5 and RTM 36 , except a few points with slightly negative values of magnitudes not exceeding 0.1 cm. We explain this by the higher smoothness of this correction compared to the harmonic correction to gravity. This results in some cancellation of positive and negative contributions to the harmonic correction generated by the mass redistribution inside the volume − as also mass redistributions in areas further away from the evaluation point still contribute.
Fifth, the harmonic correction to height anomaly is not negligible as also indicated by the histogram of Fig. 5. Peak values are 7.8 cm (RTM 5 ) and 11.5 cm (RTM 36 ). About 12% (RTM 5 ) and 23% (RTM 36 ) of all harmonic corrections are larger than 1 cm, though the majority of 77% (RTM 5 ) and 65% (RTM 36 ) of the corrections are smaller than 0.5 cm. The 95% percentiles are 1.8 cm (RTM 5 ) and 4.4 cm (RTM 36 ). Whereas the former is comparable to the reported quality of the Baltic Region and Nordic Area (NKG2015) gravimetric quasi-geoid model of 1.5-1.8 cm (95% confidence level) based on a comparison with GPS-levelling (cf. Eshagh and Berntsson 2019), the latter is significantly larger than the 95% confidence level.
Sixth, the approximation error of the harmonic correction to gravity in Forsberg and Tscherning (1981) has a non-zero mean: −7.6 mGal (RTM 5 ) and −10.2 mGal (RTM 36 ). This introduces a bias in the reduced gravity anomalies when using the harmonic correction to gravity in Forsberg and Tschern- ing ( 1981). This bias introduces long-wavelength errors in the computed quasi-geoid model. When computing a correction surface in support of GNSS-levelling, most of the bias will be absorbed.

Auvergne test area
The DTM of the Auvergne test area is based on SRTM3 (Jarvis et al. 2008). It is given for an area of 42 • N-50 • N and 2 • W-8 • E at the nodes of a 3 × 3 grid. To reduce the computational costs, we use a data window covering an area of 44 • N-49 • N and 1 • W-7 • E. We still refer to this smaller area as the Auvergne test area (cf. Fig. 6a). Likewise, the data area was reduced accordingly comprising 133,615 of the original 244,009 gravity values. Different from the Norway test area, the RTM corrections were not computed at the nodes of an equal-angular grid but at the gravity stations. The statistics of the DTM-and RTM-surfaces as well as of the tesseroid heights are given in Table 4.

Auvergne test area: harmonic correction
In the Auvergne test area, the harmonic correction had to be computed at about 59% (RTM 5 ) and 71% (RTM 36 ) of all evaluation points. These percentages are significantly larger compared to those for the Norway test area. As already found in Sect. 4.2, the smooth RTM 36 surface increases the number of evaluation points which are located below the RTM surface, i.e., for which harmonic corrections need to be computed, in the Auvergne test area. Figure 7 shows a geographic rendition of the harmonic correction to gravity of Eq. (37) and height anomaly (based on Eq. (41)), and the approxima- Fig. 7 Auvergne test area: harmonic correction to gravity, δg harm of Eq. (37), the difference between δg harm and the harmonic correction in Forsberg and Tscherning (1981), δg harm − δg FT harm , and the harmonic correction to height anomaly, δζ harm (based on Eq. (41)). The areas where topography or sea surface is located above the RTM surface are shown in white (a) δ g harm , RTM 5 (b) δ g harm , RTM 36 tion error of the harmonic correction to gravity in Forsberg and Tscherning (1981). The corresponding summary statistics are given in Table 5.
There are a few aspects which we would like to mention when comparing the results with those for the Norway test area presented in Sect. 4.2.
First, all harmonic corrections to gravity are non-negative; for the Norway test area, we found that at some evaluation points, the harmonic correction to gravity is negative for the relatively rough RTM 5 surface, which we explained by strong variations of that surface in the neighbourhood of these points. Obviously, these extreme situations do not occur in the Auvergne test area as topographic variations in this area are not that extreme.
Second, the peak approximation error of the harmonic correction to gravity in Forsberg and Tscherning (1981) appears to be very large again with −118.5 mGal (RTM 5 ) and −109.6 mGal (RTM 36 ). However, the low 95% percentiles of the absolute approximation errors, which are 4.1 mGal (RTM 5 ) and 5.2 mGal (RTM 36 ), indicate, that there is only a small number of points where such extreme values occur.
Third, the statistics of the harmonic corrections to gravity and height anomaly are overall smaller in the Auvergne test area compared to the Norway test area. This is a consequence of the moderate topography in the Auvergne test area, which implies smaller tesseroid heights compared to the Norway test area. The harmonic corrections to height anomaly are very small, though the peak values of 10.9 cm (RTM 5 ) and 13.6 cm (RTM 36 ) are even larger than the peak  (37)) and height anomaly (δζ harm , based on Eq. (41)) and of the harmonic correction to gravity in Forsberg and Tscherning (1981)  values found in the Norway test area. The large peak values may be caused by the location of the gravity stations, which in the mountains often follow the roads; correspondingly, tesseroid heights at these stations may become pretty large (remember that in the Norway test area, the evaluation points are located on a equal-angular grid and are not related to gravity stations). However, larger harmonic corrections to height anomaly only occur at a very small number of gravity stations. This is supported by the histograms of Fig. 8. About 98% (RTM 5 ) and 90% (RTM 36 ) of all corrections are smaller than 0.5 cm and only 1% (RTM 5 ) and 6% (RTM 36 ) of all corrections are larger than 1 cm. The 95% percentiles are just 0.2 cm (RTM 5 ) and 1.1 cm (RTM 36 ) compared to 1.8 cm (RTM 5 ) and 4.4 cm (RTM 36 ) in the Norway test area. From this we conclude that the harmonic correction to height anomaly is not critical when computing a quasigeoid model in the Auvergne test area. When looking at the published results of the Auvergne quasi-geoid test, standard deviations of the computed height anomalies based on a comparison with GPS-levelling data between 2.9 and 6.7 cm are reported; the median is 3.5 cm (Forsberg 2010). Only about 0.2% (RTM 5 ) and 0.8% (RTM 36 ) of the harmonic corrections for height anomaly are larger than this median. Fourth, the harmonic correction to gravity in Forsberg and Tscherning (1981) introduces a bias in the reduced gravity anomalies as already observed in the Norway test area (cf. Sect. 4.2). However, in the Auvergne test area, the magnitude of the bias is much smaller: −0.87 mGal (RTM 5 ) and −0.97 mGal (RTM 36 ). Still, it will introduce a long-wavelength bias in the computed quasi-geoid model. A suitable chosen correction surface estimated from the differ-ences between gravimetric and geometric height anomalies at a set of GNSS-levelling points will absorb most of the long-wavelength errors.

Summary and concluding remarks
We derived new (approximate) equations for the harmonic correction to gravity and potential which apply to points located below the RTM surface, and provide harmonically downward continued gravity and potential values referring to a RTM-reduced earth bounded by the RTM surface. We showed that the approximation error of the simple harmonic correction to gravity as suggested in Forsberg and Tscherning (1981) can be as large as the harmonic correction itself.
The new equations allow for the first time to compute the harmonic correction to potential and height anomaly. This is important if these types of gravity field functionals are used as data in local quasi-geoid modelling (e.g., Farahani et al. 2017;Slobbe et al. 2019). We showed that in both test areas, the harmonic correction to height anomaly can be larger than the target accuracy of today's quasi-geoid models (i.e., one centimetre) with peak values close to or even exceeding the one decimetre level at some points. Overall, however, the results indicate that the harmonic correction to height anomaly may be critical mainly in areas with strong topographic variations like the Norway test area, whereas in flat or hilly areas such as large parts of the Auvergne test area, the one-centimetre threshold will only be exceeded at a few points.
We also showed that the choice of the RTM surface has a significant impact on the magnitude of the har- monic corrections. Overall, a smooth RTM surface provides larger harmonic corrections to gravity, potential, and height anomaly. Likewise, the errors of the harmonic correction to gravity in Forsberg and Tscherning (1981) increase with increasing smoothness of the RTM surface. Therefore, we do not recommend to use very smooth RTM surfaces. From a theoretical point of view there is no need to relate the resolution of the RTM surface (i.e., its smoothness) to the spatial resolution of the GGM used in data reduction. It is sufficient to choose the smoothness of the RTM surface such that the high-frequency signals in the gravity datasets that cannot be resolved for a given data distribution and accuracy are reduced as much as possible. For some areas of interest, this may allow the choice of rougher RTM surfaces. One of the main concerns related to the harmonic correction to gravity in Forsberg and Tscherning (1981) is the bias this correction introduces in the reduced gravity anomalies. This bias causes long-wavelength errors in the computed quasi-geoid model.
The new equations for the harmonic correction to gravity and potential involve Taylor series expansions, and the dominant term of the remainder is of the order of O(ε 2 ) and O(ε 3 ), respectively. Using the results provided in Tables 2, 4, 6, and 7, we can estimate this term. For instance, the largest tesseroid height in the Norway test area is h = 1825 m (cf. Table 2, RTM 36 ), hence ε = 2.9·10 −4 . When we combine this value with the largest value of δg − RTM (Q), which is 123.3 mGal (cf. Table 6), the dominant term of the remainder of Eq. (26) is 3.1 · 10 −5 mGal. Likewise, the largest value of δV − RTM (Q) is 28.76 m 2 /s 2 (cf. Table 7, RTM 36 ); hence, the dominant term of the remainder of Eq. (32) is 6.8 · 10 −10 m 2 /s 2 . A comparison of these error estimates with the differences between the new harmonic corrections and the ones in Forsberg and Tscherning (1981) indicates how substantial the improvements are which the new harmonic corrections provide.
The complete RTM correction presented in this paper can easily be implemented in existing RTM software. The numerical complexity is about a factor of two higher than the traditional RTM correction in Forsberg and Tscherning (1981). In fact, the modified code has to be run twice; first, to compute the effect of removing the masses in the volume + ; second, to compute the effect of adding masses to the volume − . Finally, we would like to mention that the harmonic corrections to gravity, potential, and height anomaly must not be restored in the remove-compute-restore approach. They are only needed to transform values inside the masses into harmonically downward continued values and do not reflect mass re-distributions.

Data availability
The Auvergne dataset can be requested from the IAG International Service for the Geoid, https://www.isgeoid.polimi.it/. The digital elevation models for the Norway test area can be downloaded from https://eurogeographics.org/maps-for-europe/ (land topography) and https://www.gebco.net/ (bathymetry).

Open Access This article is licensed under a Creative Commons
Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Appendix A: Complete RTM correction: threestep procedure In Sect. 3, we derived the equations for the complete RTM correction in four steps. Here, we suggest a three-step procedure and show analytically that from a theoretical point of view both provide the same results. The difference with respect to the four-step procedure is that points P ∈ S which are located below the RTM surface are upward continued (in the Earth's gravity field) to the RTM surface before the masses in + are removed and masses are added to − .

Three-step procedure
Step 1 Upward-continue data at points P ∈ S located below the RTM surface to the corresponding points Q ∈ s using the free-air vertical gravity gradient; all other data are left unchanged.
where the inner gradients are defined as and P is a point on the ellipsoidal normal through P with ellipsoidal height h P + δ.
Step 2 Move the masses inside + to infinity assuming that their mass density is ρ 0 , and fill the volume − with mass to achieve a mass density of ρ 0 . Compute the effect this has on gravity and potential, respectively, at Q: Step 3 Downward-continue harmonically the data that have been upward continued in step 1 to their original locations on the Earth's surface (terrain or sea surface), i.e., from Q to P along the ellipsoidal normal through P: Using the partial derivative of Eq. (51) with respect to h, we re-write Eq. (55) as This is identical to Eq. (26), the result of the four-step procedure derived in Sect. 3.1. Inserting Eqs. (47) and (52) into Eq. (54), and observing Eqs. (46) and (51), we find Using the partial derivative of Eq. (52) with respect to h, we can simplify Eq. (57) to This is identical to Eq. (32), the result of the four-step procedure derived in Sect. 3.1. From a practical point of view, there is a minor difference between the four-step and the three-step procedure. This difference is in the upward continuation operation. Take as an example a measured gravity value. In the four-step procedure this gravity value has to be upward continued to the RTM surface using the free-air gravity gradient referring to a modified Earth which is mass-free inside + (step 2). In the three-step procedure, the upward continuation uses the freeair gravity gradient of the real Earth. Note that the gradient to be used in the upward continuation of the three-step procedure is likely rougher than the gravity gradient used in step 2 of the four-step procedure. The reason is that the removal of the masses inside + likely has a smoothing effect on the gravity gradient.

Appendix B: Semi-analytical experiments
In the first semi-analytical example, we consider the Earth to be a non-rotating homogeneous sphere of radius R and constant density ρ. The gravity measurement point P is located on the surface of that sphere, i.e., r P = R. The RTM-surface is assumed to be a sphere of radius r Q = R + h as depicted in Fig. 9a. The complete RTM correction to gravity follows from applying Eq. (26) and neglecting terms of the order O(ε 2 ): δg corr (P) = g P − g red (P) For a spherical shell, a power series expansion yields hence, Inserting the last two expressions into Eq. (59) yields With measured gravity at P, Gravity g red (P) corresponds to the harmonically downward continued gravity of our RTM-reduced earth, which is a homogenous sphere of radius R + h. Therefore, g red (P) can easily be computed exactly. For this, we consider the difference between g(P) and g red (P), which is caused by the mass M − = ρ − of a spherical shell with inner radius R and outer radius R + h, i.e., Hence, If M denotes the mass of a homogenous sphere of radius R and density ρ, the exact value of g red (P) is Equations (63) and (66) are identical up to terms O(ε 2 ). This is less than 1 μGal for any shell of thickness |h| < 9 km. The application of the derived complete RTM reduction leads to a corrected gravity value at P, which is now a boundary value of the external gravity potential caused by the mass M + = M + M − = 4 3 πρ(R + h) 3 . It is no longer an inner potential gradient and can be upward continued by the use of a free-air reduction. If necessary to apply, the free-air gradient is now based on M + . For the gravity potential of the RTMreduced earth, we find The correct analytical result of a homogeneous sphere with radius R + h and density ρ is Equations (67) and (68) agree up to terms of O(ε 4 ).
The second semi-analytical example uses a spherical cap with opening angle ψ c to model the RTM-masses in the region − (see Fig. 9b). In this more specific example, the presented precise harmonic correction formulas are also confirmed analytically.
We assume that the Earth consists of a non-rotating homogeneous sphere of radius R + h with constant density ρ with a conical mass element with symmetry axis through P and opening angle ψ c ≤ π removed (cf. Fig. 9b). The RTM-reduced earth, which results after the complete RTM correction has been applied, is a homogeneous sphere of radius R + h. With this assumption it follows directly that r P = R, i.e., the point P is located on the sphere of radius R at a depths h below the RTM surface. For this geometry, it is Hence, Eq. (35) reads g red (P) = g P − δg corr (P) The analytical expressions for the gravity effect of a spherical cap, δg sc Q , are given in Kadlec (2011, Eq. (2.106)) or Heck and Seitz (Eq. (54), 2007). Using the approximations with the Euclidean distances 1 = r 2 + r 2 1 − 2rr 1 cos ψ, 2 = r 2 + r 2 2 − 2rr 2 cos ψ , and the spherical distance ψ between the geocentric vectors of P and Q, For the last term of Eq. (73), we find ( Kadlec 2011, Eq. (2.134)) up to terms of O(ε 2 ) When inserting the last expressions into Eq. (73) the evaluations lead to the gravity value including the harmonic correction which is in accordance with the expected value (cf. Eq. (66)) on the approximation level of O(ε 2 ). For the potential the results of this thought experiment are based on Kadlec (2011, Eq. (2.98)) or Heck and Seitz (2007, Eq. (52)): (i) The gravity potential W (P) at a point P ∈ S: The expected value of the harmonically continued potential is (iii) Gravity potential at P after applying harmonic RTM reduction according Eq. (39): (iv) Items (ii) and (iii) are equal on the approximation level O(ε 3 ). The application of the derived harmonic RTM reduction leads to a corrected gravity potential at P, which is now a boundary value of the external gravity potential caused by the masses M * .

Appendix C: Complete RTM correction for the Norway test area
Figures 10 and 11 show a geographic rendition of the contributors to the complete RTM correction to gravity and potential, respectively, in the Norway test area. The statistics of the contributors are presented in Table 6 for gravity and Table 7 for potential. All contributors have distributions with a long, one-sided tail. Therefore, Tables 6 and 7 show next to range, the 25%, 75% and 95% percentiles and the median instead of mean, RMS and standard deviation. Both tables reveal that the complete RTM correction is dominated by the mass re-distribution related term, δg + RTM (P) − δg − RTM (Q) and δV + RTM (P) − δV − RTM (Q), respectively. Their contribution to the complete RTM correction is about 75% for gravity and about 90% for potential. Among the two remaining terms of the complete RTM correction to potential, which are only relevant at evaluation points located below the RTM surface, the term    The statistics were computed for all evaluation points which are located below the RTM surface. Units in mGal  Figures 12 and 13 show a geographic rendition of the contributors to the complete RTM correction to gravity and potential, respectively, over the Auvergne test area. The associated statistics are provided in Tables 8 and 9, respectively. The share of the mass redistribution related contributors to the complete RTM correction is comparable to the one found for the Norway test area, though in absolute terms, the contributors are a bit smaller now. When comparing the results for the two RTM surfaces, we notice again that the complete RTM correction and its contributors are larger for the smoother RTM 36 surface compared to the rougher RTM 5 surface for both gravity and potential.   Units are in m 2 s −2