Exact closed-form expressions for the complete RTM correction

We present exact, closed-form expressions for the complete RTM correction and the harmonic correction to disturbing potential, gravity disturbance, gravity anomaly, and height anomaly. They need to be applied in quasi-geoid modelling whenever data points are buried inside the masses after residual terrain model (RTM) reduction and analytically downward-continued functionals of the disturbing potential at the original locations of the data points are required. Compared to recent work of the authors published in this journal, no Taylor series enter the expressions and numerical instabilities of the harmonic downward continuation from the RTM surface to the Earth’s surface are avoided as are inaccuracies in the free-air upward continuation from the Earth’s surface to the RTM surface caused by a lack of precise information about higher-order derivatives of the disturbing potential. The new expressions can easily be implemented in any existing RTM software package and do not require additional computational resources. For two test areas located in western Norway and the Auvergne in France, we compute the complete RTM correction and the harmonic correction to the afore-mentioned functionals of the disturbing potential. Overall, all harmonic corrections are non-negative with maximum values of 1.54 m2/s2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{m}^2/\textrm{s}^2$$\end{document}, 263.0 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document}Gal, 263.9 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document}Gal, and 15.7 m (Norway) and 1.55 m2/s2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{m}^2/\textrm{s}^2$$\end{document}, 263.3 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document}Gal, 263.3 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document}Gal, and 15.8 cm (Auvergne) for disturbing potential, gravity disturbance, gravity anomaly, and height anomaly, respectively. The medians are 0.02 m2/s2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{m}^2/\textrm{s}^2$$\end{document}, 33.6 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document}Gal, 33.7 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document}Gal, and 0.3 cm (Norway) and 0.01 m2/s2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{m}^2/\textrm{s}^2$$\end{document}, 19.2 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document}Gal, 19.2 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document}Gal, and 0.1 cm (Auvergne). We also show that the nth Taylor polynomials used in the recent work of the authors published in this journal may have large remainders depending on the topography in the vicinity of the evaluation point no matter how n is chosen. Finally, we show that the commonly used expression for the harmonic correction to gravity anomaly introduced in 1981 is almost exact, though it was derived along a completely different line of reasoning. The errors do not exceed 49 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document}Gal in both test areas. Moreover, the errors have a negligible impact on the computed height anomalies in one-centimetre quasi-geoid modelling, as the mean error does not exceed a few μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document}Gal in both test areas.


Introduction
The residual terrain model (RTM) reduction was introduced in Forsberg and Tscherning (1981) as a method to deal with topography in gravity field modelling. Basically, the RTM reduction aims at removing from the data the contribution of short-wavelengths topographic signals relative with respect to a smooth mean elevation surface, here referred to as the RTM surface. In this way, unwanted high-frequency topographic signals present in the data are largely removed depending among others on the quality and resolution of the topographic information (i.e. digital elevation model and/or bathymetry model, mass density model). After application of the RTM reduction, stations above the RTM surface are essentially in free-air while stations below the RTM surface are buried inside the masses. The latter implies that the RTMreduced geopotential functionals at stations located below the RTM surface represent values inside the masses. Before these RTM-reduced internal geopotential functionals are used in quasi-geoid modelling, Forsberg and Tscherning (1981) suggested to apply a correction, referred to as the harmonic correction, which provides analytically downward-continued functionals at the original location of the data points.
The expression for the harmonic correction to gravity anomaly suggested in Forsberg and Tscherning (1981), though very simple to evaluate, was known to be an approximation, which at that time was sufficient as no highresolution, precise digital elevation models were available. Forsberg and Tscherning (1981) did not derive expressions for the harmonic correction to other functionals of the disturbing potential. However, as far as disturbing potential and height anomaly are concerned, they argued that the corresponding harmonic corrections are likely very small and may be neglected in most cases. Since then, many attempts were made to address the problem of the harmonic correction (cf. Klees et al. 2022 for a review).
Recently, Klees et al. (2022) provided for the first time expressions for the harmonic correction to disturbing potential and height anomaly and improved expressions for the radial derivative of the disturbing potential. Using topographic data over a part of western Norway and the Auvergne in France, they computed and analysed the harmonic correction to these functionals. In fact, their expressions for the harmonic correction are also approximations, as they were based on nth Taylor polynomials. Based on an estimated upper bound of the remainders of the nth Taylor polynomials, ibid. suggested the choice n = 2 for disturbing potential and height anomaly and n = 1 for the radial derivative of the disturbing potential. Using the new expressions, they showed that in areas of strong topographic variations (e.g. around the Norwegian fjords), the differences with respect to the harmonic correction to gravity anomaly in Forsberg and Tscherning (1981) can be as large as 100 mGal. Moreover, as these differences do not average out over the data area, they may introduce significant long-wavelength errors in the computed quasi-geoid model. They also showed that the harmonic correction to disturbing potential and height anomaly only need to be applied in areas of strong topographic variations with peak values of 1 m 2 /s 2 (disturbing potential) and 0.1 m (height anomaly) in the Norway and Auvergne test areas, whereas in the flat and hilly parts of these areas, they are mostly below 0.1 m 2 /s 2 (disturbing potential) and 0.01 m (height anomaly).
The approach suggested in Klees et al. (2022) allows to compute the complete RTM correction (i.e. the sum of RTM reduction and harmonic correction) to any linear functional of the disturbing potential. It involves four steps (cf. Fig. 1): 1. Masses inside the volume + are removed, and the functional is corrected for the effect this has.
2. The functional obtained in step 1 is upward-continued in free air from its original location on the Earth's surface to the RTM surface along the ellipsoidal normal. 3. Masses are added to the volume − to obtain a constant mass density in − , and the functional from step 2 is corrected accordingly. 4. The functional obtained in step 3 is downward-continued analytically from the RTM surface to the original location on the Earth's surface along the ellipsoidal normal providing an analytically downward-continued functional at the original location. This functional is ready to be used in quasi-geoid modelling. Klees et al. (2022) used truncated Taylor series for the freeair upward continuation of step 2 and the analytic downward continuation of step 4, respectively. For each remainder, they provided upper bounds based on the exterior gravity field of a homogeneous mean Earth sphere. However, we will show in Appendix B using a simple example that the upper bounds provided in Klees et al. (2022) are not really upper bounds everywhere. In particular, in deep, narrow valleys, they may underestimate the magnitude of the remainders significantly. Consequently, some results and conclusions presented in Klees et al. (2022) became wrong and need to be corrected. Moreover, we also show in Appendix B that this problem is not solved by just increasing the order n of the Taylor polynomials. Therefore, we need another approach for step 2 and step 4, which does not involve Taylor series.
In this paper, we will provide exact, closed-form expressions for the complete RTM correction to disturbing potential, gravity disturbance, gravity anomaly, and height anomaly. They are easy to implement in any existing RTM software package and do not require additional computational resources. In Sect. 2, we will derive the exact, closed-form expression for disturbing potential. A key quantity arising in this expression is the harmonic downward continuation of the gravitational field of the masses added to volume − (cf. Fig. 1) to the location of a data point, which, after step 3, is buried in the masses. An exact, closed-form expression for this harmonic downward continuation is derived in Sect. 3. The results of Sects. 2 and 3 are combined in Sect. 4 providing the final exact, closed-form expression for the complete RTM correction to disturbing potential and also its radial derivative. Likewise, exact, closed-form expressions are provided for gravity disturbance, gravity anomaly, and height anomaly using their relation with the disturbing potential. In Sect. 5, the new expressions will be applied to the two test areas already considered in Klees et al. (2022). Moreover, we will identify the results and conclusions in Klees et al. (2022) which are wrong due to the underestimation of the magnitude of the remainders of the nth Taylor polynomials used in step 2 and step 4 of the complete RTM correction approach. In Sect. 6, we close with a summary and some final remarks. To derive the new expression for the complete RTM correction to disturbing potential, we start with the four-step procedure suggested in Klees et al. (2022), but use a more compact notation for a better readability. Let T be the disturbing potential, P a point on the Earth's surface, and Q the intersection of the ellipsoidal normal through P and the RTM surface (cf. Fig. 1). If P is located above the RTM surface, the complete RTM correction reduces to the RTM reduction introduced in Forsberg and Tscherning (1981). After RTM reduction, the point P is hanging in free-air inside the harmonic domain of the RTMreduced disturbing potential. Here, we are only interested in points P, which are located below the RTM surface. These points are buried in the masses, which were added to the volume − in step 3. We denote by V + the gravitational potential of the masses in + and by V − the gravitational potential of the masses added to volume − (cf. Fig. 1). Note that if our data area does not contain water surfaces, − is mass-free before step 3 (cf. Fig. 2 and Table 1 in Klees et al. (2022)). Then, the masses added to − in step 3 are identical Step 1: Move the masses in + to infinity and correct T P for the gravitational effect this has: Step 2: Upward continue in free-air T + P of Eq.
(1) to the point Q on the RTM surface along the ellipsoidal normal through P: Step 3: Add mass to the volume − so that the mass density equals some pre-selected value and correct T + Q of Eq. (2) for the gravitational effect this has: Step 4: Downward-continue harmonically T red from point Q on the RTM surface to point P on the Earth's surface along the ellipsoidal normal through P: Note that the upper index "UC" in Eq.
(2) and "HDC" in Eq. (4) means "free-air upward continuation" and "harmonic downward continuation", respectively. When combining Eqs. (1)-(4), we can write T red,HDC P as the difference between T P and the complete RTM correction, denoted (T ) cRTM We find it important to mention that due to step 4, T red,HDC P is different from T red P , the value of T red at P inside the masses. Remember that after step 3, T red is actually harmonic outside the RTM surface ("actually" means we ignore that the DEM is just an approximation to the Earth's surface and the density of the masses in + and of the masses in − (if there are any), are not exactly known). In step 4, T red is harmonically downward-continued from point Q on the RTM surface through the masses to point P on the Earth's surface, which gives T red,HDC P . Using the definition of the RTM reduction according to Forsberg and Tscherning (1981), Equation (5) is re-written as The term in brackets on the right-hand side of Eq. (7) is the exact expression for the harmonic correction, denoted (T ) harm P . Hence, we may write with V − P , which does not appear in the expression for the complete RTM correction appears by construction in the RTM reduction in Eq. (6) and the harmonic correction of Eq. (9). Note that V − P is the gravitational potential of the masses added to volume − in step 3, evaluated at P. Though P is located inside the masses after step 3, P is also a point located on the boundary of the harmonic domain of the gravitational potential of the masses added to volume − .
In Klees et al. (2022), the effect of the free-air upward continuation (UC) from point P to point Q in step 2, δT UC P→Q , and the effect of the harmonic downward continuation (HDC) from point Q to point P in step 4, δT HDC Q→P , were considered separately. For each effect, a nth Taylor polynomial about P (δT UC P→Q ) and Q (δT HDC Q→P ) was used. In Appendix B, we have chosen a simple mass configuration to demonstrate that the remainder of each nth Taylor polynomial used in Klees et al. (2022) may introduce errors in the computed complete RTM correction, which are much larger than the upper bounds in Klees et al. (2022) suggest and also much larger than the measurement error. This applies in particular to the effect of the harmonic downward continuation, δT HDC Q→P , at points P located in deep, narrow valleys, i.e. at locations where the complete RTM corrections attain large magnitudes. The remainder of the harmonic downward continuation in step 4 does not decrease when increasing n. Moreover, in deep, narrow valleys, the upward continuation in step 2 may also be prone of errors, because derivatives of the disturbing potential may differ significantly from their nominal values.
We will show that Taylor series for the upward continuation of step 2 and the harmonic downward continuation of step 4 can be avoided completely. To do so, we exploit the superposition principle of gravitating masses, which allows the combination of step 2 and step 4 into a single step, and provides closed-form expressions for the complete RTM correction. Key are three observations. First, the upward continuation of step 2 is done in free-air in the gravitational field T + (cf. Eq. (2)). Second, the point P is located on the boundary of the harmonic domain of T + as T + is the disturbing potential before masses are added to volume − . Third, the disturbing potential after step 3, i.e. after masses has been added to volume − , T red = T + + V − , is harmonic outside the RTM surface. Hence, the continuity of T + across the boundary of the harmonic domain allows us to write T red,HDC where V −,HDC P is the harmonic downward continuation of the gravitational potential V − from its upper harmonic domain through the masses to point P. Using Eq. (10), the effect of the harmonic downward continuation of step 4, δT HDC Q→P , can be written as In Eq. (11), we move δT UC P→Q to the left-hand side, which provides an expression for the joint effect of the free-air upward continuation of step 2 and the harmonic downward continuation of step 4: Equation (12) is the first main result of this paper. The sum appears in both the expression for the complete RTM correction and the harmonic correction. Therefore, with Eq. (12) we may write for the complete RTM correction, (T ) cRTM P , Eq. (8), and the harmonic correction, (T ) harm P , Eq. (9) at points P located below the RTM surface: Equation (13) consists of two terms. The term V + P is the gravitational potential at P of the masses in + . The term V −,HDC P is the harmonic downward continuation to point P of the exterior gravitational potential of the masses added to volume − , where "exterior" refers to the region above the RTM surface.
The harmonic correction to disturbing potential of Eq. (14) also comprises two terms, one of them is V −,HDC P , which already appeared in Eq. (13). The other one, V − P , is the gravitational potential at P of the masses added to volume − . Note that P is located on the lower boundary of the harmonic domain of V − . V + P and V − P are computed when computing the RTM reduction (cf. Eq. (6)). To do so, the volumes + and − are discretized using rectangular prisms and/or tesseroids. For prisms, a closed-form expression exists (cf. Nagy et al. 2001Nagy et al. , 2002. For tesseroids, a truncated Taylor series is used (cf. Heck and Seitz 2007;Grombein et al. 2013).
3 The computation of V −,HDC P V −,HDC P plays a key role in the complete RTM correction and the harmonic correction according to Eqs. (13) and (14). In principle, it could be computed using an nth Taylor polynomial about Q. However, as shown in Appendix B, this does not provide an accurate enough approximation to V −,HDC P , in particular not in deep, narrow valleys, no matter how n is chosen. Finding a suitable, preferably closed-form expression for V −,HDC P is not straightforward and requires some additional considerations.
Suppose we cut out from − a cylinder of radius R and height h = h Q − h P , and a symmetry axis identical to the line Q − P (cf. Fig. 2). The volume cut out in this way is denoted − h , the remaining volume is denoted δ − , and the volume of the cylinder is denoted c . Then, The difference ε := − h − c is the volume between i) the top of the cylinder and the RTM surface, and ii) the bottom of the cylinder and the Earth's surface (cf. Fig 2). We find it important to mention here that we can make this difference arbitrarily small by choosing a small enough radius R. Now, according to the superposition principle of gravitating masses, we can write V −,HDC P as the sum of the contribution of the masses in δ − and − h , i.e.
Because P is located on the boundary of the harmonic domain of the masses in δ − , the first term on the right-hand side of Eq. (16) is Note that V − P − h is the gravitational potential at point P of the masses which were cut out. Using Eqs. (41) and (43) in Appendix A, it is Inserting Eqs. (17)-(19) in Eq. (16), we find However, according to Eqs. (41) and (43) independent of R. Hence, The size of the term {·} on the right-hand side of Eq. (22) is controlled by the size of R. In particular, Hence, we can make this term arbitrarily small by choosing a sufficiently small R. At the same time, the right-hand side of Eq. (21) does not depend on R. Hence, when choosing a sufficiently small R, Eq. (22) is actually equal to and can be considered as being exact. Equation (24) is the second main result of this paper.

Exact, closed-form expressions
Using Eq. (24), we obtain the following expression for the complete RTM correction to disturbing potential, Eq. (13), and the harmonic correction to disturbing potential, Eq. (14): Both equations are of closed-form and, as explained in Sect. 3, actually exact. Following the same line of reasoning, we can easily find exact, closed-form expressions for the complete RTM correction and harmonic correction to the radial derivative of the disturbing potential. With (e.g. Kadlec 2011 with the remark that in Eq. 2.52 of ibid. it must read −h 2 instead of h 2 ), which is also independent of R, we find To find similar expressions for the complete RTM correction and the harmonic correction to gravity disturbance, δg, gravity anomaly, g, and height anomaly, ζ is straightforward. First, we remember the relationship between these functionals and the disturbing potential (we use the spherical approximation for convenience): where r is the spherical distance to and γ the normal gravity at the telluroid point associated with the evaluation point. Second, we define Then, using Eqs. (25)-(28), we find: We find it important to mention that in each of the Eqs. (35), (37), and (39) the term in brackets on the right-hand side is equal to the corresponding RTM reduction. In practice, only the expressions for the complete RTM correction are relevant. However, the expressions for the harmonic corrections are of interest if one is interested in the error of the harmonic corrections to gravity anomaly and gravity disturbance in Forsberg and Tscherning (1981) and (Klees et al. 2022) or in the magnitude of the harmonic correction to disturbing potential and height anomaly. Remember that Eq. (38) is actually an exact expression for the harmonic correction to gravity anomaly. It differs from the harmonic correction to gravity anomaly suggested in Forsberg and Tscherning (1981) by just the term −4π Gρ h 2 r . For the Norway and Auvergne test areas (see Sect. 5), extreme tesseroid heights do not exceed 2000 m. Then, the term 4π Gρ h 2 r has a magnitude not exceeding 0.14 mGal when assuming a mass density ρ = 2670 kg/m 3 . This result is very astonishing because the harmonic correction to gravity anomaly in Forsberg and Tscherning (1981) was derived along a completely different line of reasoning.

Application to the Norway and Auvergne test areas
We apply the expressions for the complete RTM correction and the harmonic correction to the Norway and Auvergne test areas. For details about these areas including statistics of the topography, we refer to Klees et al. (2022). Whereas (Klees et al. 2022) used a low-resolution and a high-resolution RTM surface, we will present the results for the low-resolution RTM surface, which (Klees et al. 2022) called RTM 36 . RTM 36 mimics a reference gravity field complete to degree 300. It was computed as by applying a spherical Gaussian lowpass filter to the high-resolution DEM, which was truncated at its half width of 36 . We find it important to mention that for the Norway test area, the evaluation points are the nodes of an 18 × 18 grid, which formed a subgrid of the DEM, whereas for the Auvergne test area, the evaluation points are identical to the gravity stations. Figure 3 shows a geographic rendition, Fig. 4 the histograms, and Table 1 some statistics of the complete RTM correction to gravity anomaly and height anomaly for the Norway and Auvergne test areas. Only evaluation points below the RTM surface have been included (222, 892 for the Norway test area, which corresponds to 28% of all points; 95, 360 for the Auvergne test area, which corresponds to 71% of all points). The corrections have large peak magnitudes of 377.2 mGal and 1.159 m (Norway test area) and 163.5 mGal and 2.128 m (Auvergne test area). The different topographic regimes in both test areas explain the differences in the complete RTM correction both in terms of spatial pattern, distribution, and range. The bi-model distribution of the complete RTM correction to height anomaly in the Auvergne test area can be explained by the distinguished topography with rather flat regions left to a line from south-west to northeast and moderate to high mountain ranges right to that line (cf. Fig. 6 in Klees et al. (2022)).
We do not show corresponding information for gravity disturbances and disturbing potential. The reason why we left out the results for gravity disturbances is that the complete RTM correction does not differ much from that for gravity anomalies as the histogram of the differences in Fig. 5 reveals. In both test areas, the maximum difference is 48 µGal, and the median is below 1 µGal. The complete RTM correction to disturbing potential is not shown as it is actually a scaled complete RTM correction to height anomaly. Figures 6 and 7 show a geographic rendition and the histograms of the harmonic correction to gravity anomaly and height anomaly, respectively, for the Norway and the Auvergne test areas. Table 2 complements this information with some statistics of the harmonic corrections to disturbing potential, gravity disturbance, gravity anomaly, and height anomaly. All harmonic corrections are non-negative. Moreover, peak values of the harmonic correction to grav-ity anomaly and height anomaly are comparable in both test areas (263.9 mGal and 15.7 cm in the Norway test area versus 263.3 mGal and 15.8 cm in the Auvergne test area). However, there are significant differences in the distribution of the harmonic corrections with overall larger values in the Norway test area compared to the Auvergne test area (cf. Fig. 7). For instance, the 95% percentiles are 161 mGal and 5.9 cm in the Norway test area compared to 72.7 mGal and 1.2 cm in the Auvergne test area.
The majority (60% in the Norway test area and 89% in the Auvergne test area) of the harmonic corrections to height anomaly do not exceed 0.5 cm. Only 25% (Norway) and 6% (Auvergne) are larger than 1 cm, and only 5% are larger than 5.9 cm (Norway) and 1.2 cm (Auvergne). Values of several centimetres are only attained at a relatively low number of points. For instance, among the 222, 892 (Norway) and 95, 360 (Auvergne) evaluation points located below the RTM surface, there are only 13, 818 (Norway) and 651 (Auvergne) evaluation points with corrections exceeding 5 cm, and 2, 461 (Norway) and 188 (Auvergne) evaluation points with corrections exceeding 10 cm.
Differences between the harmonic correction to gravity anomaly and gravity disturbance are small as one can expect when comparing Eq. (38) with Eq. (36). The maximum difference is 49 µGal, and only 10% exceed 10 µGal in both test areas. These differences are at the same time identical to the approximation error of the harmonic correction to gravity anomaly in Forsberg and Tscherning (1981). The mean approximation error is just 3 µGal (Norway) and 0.9 µGal (Auvergne). Basically, a nonzero mean error of the harmonic correction to gravity anomaly in Forsberg and Tscherning (1981) may cause a bias in the computed height anomalies. The magnitude of the bias depends on the area covered by the biased gravity anomaly dataset. However, seen the small magnitude of the mean error, one can expect negligible biases if one aims at a one-centimetre quasi-geoid model. For instance, assuming that the size of the area is 500 × 500 km 2 , a bias of 3 µGal causes a maximum height anomaly error of less than 1 mm (cf. Klees and Slobbe 2018). In reality, the bias is even smaller as it only occurs at stations which are located below the RTM surface.

Rectification of the results in Klees et al. (2022)
The harmonic correction to gravity disturbance and disturbing potential published in Klees et al. (2022) were computed using nth Taylor polynomials for step 2 and step 4 of the fourstep procedure with n = 1 for gravity disturbance and n = 2 for disturbing potential. The upper bounds of the remainders of the nth Taylor polynomials provided in Klees et al. (2022) were estimated using a homogeneous spherical Earth. These  Histogram of the differences between the complete RTM correction to gravity anomaly and gravity disturbance at evaluation points located below the RTM surface upper bounds appeared to underestimate the magnitude of the remainders. Therefore, the statistics of the harmonic correction to gravity disturbance and disturbing potential in Klees et al. (2022) are not correct. This concerns in particular the extreme values and, correspondingly, the range. The largest errors are attained in deep, narrow valleys, which are numerous in the Norway test area and less frequent in the Auvergne test area. The same applies to the complete RTM correction in Klees et al. (2022), because ibid. computed it as sum of RTM reduction and harmonic correction. Here, we provide a rectification of the main results and conclusions in Klees et al. (2022): • Klees et al. (2022) reported maximum errors of the harmonic correction to gravity anomaly suggested in Forsberg and Tscherning (1981) as large as the harmonic correction itself in extreme cases. This result is wrong. In reality, the error of the harmonic correction to gravity anomaly suggested in Forsberg and Tscherning (1981)    in Forsberg and Tscherning (1981) is exact for gravity disturbances. • Klees et al. (2022) reported harmonic corrections to gravity disturbance which may be negative at some locations.
In reality, all harmonic corrections are non-negative. • Klees et al. (2022) reported large nonzero mean errors for the harmonic correction to gravity anomaly suggested in Forsberg and Tscherning (1981), in particular for the Norway test area, which may bias the computed quasi-geoid model. This result is wrong. In reality, the mean errors for the harmonic correction to gravity anomaly suggested in Forsberg and Tscherning (1981) are very small in the Norway and Auvergne test areas, not exceeding a few µGal. The corresponding bias in the computed quasigeoid model is negligible if one aims at a one-centimetre quasi-geoid model.

Summary and final remarks
We derived exact, closed-form expressions for the complete RTM correction and the harmonic correction to disturbing potential, gravity disturbance, gravity anomaly, and height anomaly using the four-step procedure suggested in Klees et al. (2022). The derivation was motivated by the fact that the expressions in Klees et al. (2022), which were based on the use of nth Taylor polynomials, appeared to be poor approximations for particular evaluation points and topographies. Particularly, the harmonic downward continuation in step 4 appeared to be unstable in some cases, which rendered the use of any nth Taylor polynomial useless, and required a Taylor-series-free approach. Key towards the new expressions was the observation that the free-air harmonic upward continuation in step 2 and the harmonic downward continuation in step 4 cancel out each other to a large extent; what is left is identical to a harmonic downward continuation in the gravitational field of the masses added to volume − . For this harmonic downward continuation, we derived exact, closed-form expressions for the disturbing potential and its radial derivative using the superposition principle of gravitating masses. This leads to surprisingly simple expressions for the complete RTM correction and the harmonic correction to disturbing potential, gravity disturbance, gravity anomaly, and height anomaly. Astonishing is the observation that the harmonic correction to gravity anomaly introduced in Forsberg and Tscherning (1981) is almost exact, with maximum errors in the two test areas not exceeding 49 µ Gal. Even more, it is exact if interpreted as the harmonic correction to gravity disturbance. The new expressions for the complete RTM correction and the harmonic correction to disturbing potential, gravity disturbance, gravity anomaly and height anomaly are easy to implement in any existing RTM software package and do not require more resources than the computation of the harmonic correction to gravity anomaly suggested in Forsberg and Tscherning (1981).
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
Here, we verify the expressions for the complete RTM correction and the harmonic correction derived in Sect. 4 using two synthetic mass distributions for which closed-form expressions for all relevant quantities of the four-step procedure in Klees et al. (2022) are known. Without loss of generality, we set the disturbing potential equal to the gravitational potential of the mass distribution and assume that the volume + is free of masses, i.e. V + = 0 and ∂ r V + = 0. First, we verify the expressions for the complete RTM correction and the harmonic correction to disturbing potential. The mass distribution is given by a large cylinder (radius R l , height h l , and constant mass density ρ) with a small cylinder (radius R s < R l , height h s < h l , and constant mass density ρ) cut out. The small cylinder represents − . The two cylinders share the same symmetry axis and are aligned along the top (see Fig. 8). In step 3 of the four-step procedure, the volume − is filled with masses, which corresponds to putting back the small cylinder. Hence, T red is identical to the gravitational potential of the large cylinder. We use well-known expressions for the gravitational potential of a homogeneous cylinder for points located on the symmetry axis, which can be found in, e.g. Heiskanen andMoritz (1967) and(Kadlec 2011). In the following, we need to know the expression for the gravitational potential of a homogeneous cylinder (radius R, height h, constant mass density ρ) at three selected points on the symmetry axis: above the cylinder, inside the cylinder, and below the cylinder. We choose the symmetry axis as the z-axis, put the origin z = 0 at the top of the cylinder, and assume that the positive z-axis points upwards (cf. Fig 9). Then (cf. Kadlec (2011)), where Moreover, due to the continuity of the gravitational potential in R 3 , it is Point P and point Q are located on the symmetry axis with z P = −h s and z Q = 0, respectively. It is (remember that in our example there are no masses in + , i.e. V + = 0 ⇒ T = T + ): T red,HDC Note that T red,HDC P is the harmonic downward continuation of T red from above the large cylinder to point P. From these equations, we can compute the complete RTM correction and the harmonic correction as Equations (58) and (59) now need to be compared with the corresponding equations derived in Sect. 4, which for the given mass configuration are The results are shown in Table 3. They confirm that the expressions for the complete RTM correction and the harmonic correction to disturbing potential of Sect. 4 are correct. Next, we verify the expressions for the complete RTM correction and harmonic correction to the radial derivative of the disturbing potential. We consider a synthetic Earth, which is a homogeneous sphere of radius R and gravitational constant G M. The gravitational potential of the synthetic Earth is set equal to the disturbing potential. The RTM surface is assumed to be the surface r = R + h, where r is the radial distance from the centre of mass of the synthetic Earth. Hence, the volume − is a spherical shell with inner radius R and outer radius R + h. In step 3 of the four-step procedure of Sect. 2, we fill the mass-free volume − with mass of density ρ. Let the total mass in − be denoted M − . There are no masses in + , i.e. ∂ r V + = 0 ⇒ ∂ r T P = ∂ r T + P . Note that r P = R and r Q = R + h. Using afore-mentioned mass distributions, the following expressions apply: Using these expressions, we can compute the complete RTM correction, (∂ r T ) cRTM P and the harmonic correction, Note that in the last equation, we used the fact that ∂ r V − P = 0. Remember that the radial derivative of the gravitational potential of a homogeneous spherical shell is continuous in R 3 and zero inside the spherical shell.
Equations (67) and (68) need to be identical with the expressions for the complete RTM correction to ∂ r T , Eq. (27), and the harmonic correction to ∂ r T , Eq. (28), which are for the given mass configuration equal to (remember that ∂ r V + P = 0) The results in Table 4 provide the evidence, that this is indeed so, i.e. the expressions for the complete RTM correction and the harmonic correction to the radial derivative of the disturbing potential of Sect. 4 are correct.

Appendix B
We want to demonstrate that using a Taylor series for the harmonic downward continuation of step 4 as suggested in Klees et al. (2022) may be unstable in some cases with higher order terms increasing rapidly with alternating signs. Moreover, we show that using a nth Taylor polynomial with n = 1 for gravity anomaly and gravity disturbance, and n = 2 for disturbing potential and height anomaly as suggested in This also implies that the upper bounds provided in Klees et al. (2022) are not correct. In particular, in deep, narrow valleys (i.e. in areas where harmonic corrections take up large values), they may underestimate the magnitude of the remainder significantly. Note that the problem is not solved automatically by computing the complete RTM correction, simply because the harmonic downward continuation in step 4 is also part of the complete RTM correction. Combining step 2 and step 4 also does not solve this problem, as an analytic downward continuation through the masses in − also appears in the sum of step 2 and step 4. The solution to this problem, which was provided in Sect. 4 indeed solves this problem.
We consider a single cylinder of radius R and height h (cf. Fig. 9). As in Appendix A, the z-axis is identical to the symmetry axis, the top of the cylinder is z = 0 and the positive z-axis points upwards. Equations (41)-(43) describe the gravitational potential of the cylinder for points on the symmetry axis above, inside and below the cylinder, respectively. We consider the harmonic downward continuation of the gravitational potential V z>0 (R, h, z) from a point Q on top of the cylinder (i.e. z Q = 0) to a point P at the bottom of the cylinder (i.e. z P = −h). The harmonic downward continuation of the gravitational potential above the top of the cylinder to the point P is given by V z>0 (R, h, z P ). Klees et al. (2022) used the nth Taylor polynomial about z = z Q = 0 to compute the harmonic downward continuation to point P as They claim that n = 2 provides a sufficiently accurate approximation to V HDC P . We will show using the simple example of a homogeneous cylinder that V HDC,n P may dif- fer significantly from V HDC P not only for n = 2, but also for larger values of n. Table 5 shows the results for various choices of R, h, and n. Obviously, V HDC,n P is a poor approximation to V HDC P if h/R is large (e.g. column 2 and 3 in Table 5). If the ratio h/R decreases, the approximation error decreases with increasing n. An analysis of the term in the Taylor series reveals that it comprises a factor h R m−1 , which dominates this term. Hence, for a cylinder, we can expect no convergence of V HDC,n P if h/R > 1. On the other hand, if h/R < 1, V HDC,n P may provide a good approximation to V HDC P already for moderate values n. For practical applications, this result implies that the harmonic downward continuation from Q on the RTM surface to P on the Earth's surface using a nth Taylor polynomial does not provide a good approximation in deep, narrow valleys (the equivalent of h/R > 1 for a cylinder). Only if the width of the valley is greater than twice its depth (corresponding to h/R < 1 for the cylinder), the nth Taylor polynomial may provide a reasonable (relative error on the order of a few percent) approximation already for low n, e.g. n = 2 as suggested in Klees et al. (2022). In general, however, a nth Taylor polynomial should not be used for the harmonic downward continuation V HDC P . The same applies to the analytic downward continuation of linear functionals of the disturbing potential such as gravity anomalies, gravity disturbances, and height anomalies.