The determination of the geoid in South Korea is a national imperative for the modernization of height datums, specifically the orthometric height and the dynamic height, that are used to monitor hydrological systems and environments with accuracy and easy revision, if necessary. The geometric heights above a reference ellipsoid, routinely obtained by GPS, lead immediately to vertical control with respect to the geoid for hydrological purposes if the geoid height above the ellipsoid is known accurately. The geoid height is determined from gravimetric data, traditionally ground data, but in recent times also from airborne data. This paper illustrates the basic concepts for combining these two types of data and gives a preliminary performance assessment of either set or their combination for the determination of the geoid in South Korea. It is shown that the most critical aspect of the combination is the gravitational effect of the topographic masses above the geoid, which, if not properly taken into account, introduces a significant bias of about 8 mgal in the gravity anomalies, and which can lead to geoid height bias errors of up to 10 cm. It is further confirmed and concluded that achieving better than 5 cm precision in geoid heights from gravimetry remains a challenge that can be surmounted only with the proper combination of terrestrial and airborne data, thus realizing higher data resolution over most of South Korea than currently available solely from the airborne data. The accurate determination of the geoid, or level reference surface for heights, is a major goal of infrastructure modernization in many countries. For example, knowing the geoid height, N , above a reference ellipsoid everywhere in a country permits easy determination of orthometric heights, H , using the Global Positioning System (GPS), or any other Global Navigation System of Satellites (GNSS), thus avoiding labor intensive spirit leveling procedures that were traditionally used to establish vertical control: where h is the GPS-measured height above the ellipsoid (Hofmann-Wellenhof and Moritz, 2005). With an accurate model for N , heights relative to the geoid may be determined with uniform accuracy anywhere a GPS receiver antenna can be placed; and, they may be easily and quickly reestablished in areas subject to surface deformation, erosion, or any sudden change caused by natural disasters such as flooding or earthquakes. Determining the geoid heights requires an accurate, detailed, and spatially extensive data base of gravity values. Traditional measurements of gravity were and still are obtained with instruments placed at surface stations; in this paper these are called terrestrial data. In easily accessible areas such stations may be fairly regularly distributed. However, often in topographically rugged regions, they are spaced along roads, which also typically follow valleys; and, in these cases they are naturally sparsely distributed. Over the last few decades, airborne gravimetry has matured to the point that regional surveys for geoid determination are readily accomplished with uniform accuracy and spatial resolution (Forsberg and Olesen, 2010). These data are somewhat less accurate than terrestrial measurements, but cover otherwise inaccessible areas such as mountainous terrain and coastal regions, where ship surveys also may not be possible or are very expensive. Moreover, in these areas airborne gravimetry fills the resolution gap that is left by the satellite gravity missions, such as GOCE (Gravity Field and Steady-State Ocean Circulation Explorer, Rummel et al ., 2011), which achieve a maximum resolution of about 85 km. While airborne gravimetry may be an economical methodology to build a data base for geoid determination, existing and new terrestrial measurements should be included, not just for data base enhancement, but also to help calibrate and validate airborne data. Certainly both types of data must refer to the same gravitational field and thus must be consistent at the level of the airborne measurement accuracy, about 1-2 mgal (1 mgal =10 ⁻⁵ m/s ² ). However, experience shows that initially, at least, this is not the case. For the recent airborne survey completed over South Korea in 2008-2009 (Bae et al ., 2012), the differences between terrestrial and airborne gravity data at common geographic locations throughout the country have a mean value of about 8 mgal and a standard deviation of about 12 mgal (see Section 4, below). This observed bias between the terrestrial and airborne data in South Korea is not due to the attenuation of the gravity field with altitude. Indeed, it is just the opposite: the airborne data are significantly greater in value than the terrestrial data in many cases. Leaving a significant bias in the gravity data leads to a bias in the geoid determination, and resolving the bias is, therefore, an important part in the combination of the two data types. It is found that the gravitational effect of the topographic masses above the geoid is largely responsible for the bias. However, reducing the data for this gravitational effect requires special care, since the airborne and terrestrial data have fundamentally different characteristics. That is, aside from the obvious difference in vertical height of the measurement surface, the data also have different spatial distribution and spectral content. In this paper, we briefly review the geoid height determination from gravimetric data, discuss the basic concepts of computing the terrain effect, and on the basis of a quantitative assessment propose methods to combine the terrestrial and airborne data. This appraisal includes a preliminary analysis of the subsequent geoid height errors, which also confirms that the available combined gravimetric data have not yet reached the resolution and accuracy needed to obtain better than 5 cm accuracy in the geoid height. The detailed, precise geoid height is determined principally from a distribution of gravity data on or near the Earth’s surface. The long-wavelength geoid height, on the other hand, typically is based on a high-degree reference model, determined either from a combination of satellite orbit analysis and global gravity data or, more recently, from in situ satellite measurements of the gravitational field. This long-wavelength model is required since theory otherwise demands that a global distribution of gravity data enter into the computations. In addition, the local topography can be used partly to aid in the estimation of the shorter wavelengths of the geoid height; and, its gravitational effect also must be treated separately in order to fulfill the theoretical condition that the data surface is a level surface (e.g., the geoid, itself). In summary, the equation in spherical approximation that yields the geoid height, N , from gravity anomalies, Δg , is (Hofmann- Wellenhof and Moritz, 2005; Wang et al ., 2012; Huang and Véronneau, 2013) where N ₀ is a constant offset of the local geoid from the global geoid, γ ₀ is an average value of gravity; R is the mean radius of the geoid approximated as a sphere; the point of computation is ξ = ( θ ,λ ), represented by spherical polar coordinates; σ is the unit sphere with dσ = sin θ' dθ' dλ' ; S ( ψ ) is Stokes’s function with ψ the central angle between the computation and integration points; the reference gravity anomaly and geoid height are spherical harmonic models, respectively. where GM is Newton’s gravitational constant times Earth’s total mass; a is the radius to which the harmonic coefficients, C _nm , of a reference field refer; ( ξ ) is the fully normalized spherical harmonic function; and, δN _TE is the indirect effect associated with the processing of the terrain effect, δg _TE . The gravity anomalies, Δg , are differences between measured gravity and a normal gravity based on a reference, or normal potential, and are known as free-air anomalies. With the inclusion of the reference field, equations (3) and (4) , the integral in equation (2) may be truncated to a region, like South Korea. This also offers the ability to modify Stokes’s function in order to further reduce the error committed in the truncation of the integral (Jekeli, 1981; Featherstone, 2013). The gravity anomalies refer to the Earth’s surface, or to the geoid with a free-air reduction, where the latter usually involves a simple linear approximation to the vertical gradient of gravity. As such they represent the needed boundary values that justify the integral solution in equation (2). However, this solution further requires that the boundary actually bounds all masses. To achieve this condition, the terrain is mathematically moved inside (or onto) the boundary with a corresponding adjustment to the data values. This adjustment is known as the direct terrain effect. It comprises the computation of the Bouguer anomaly, which is the free-air anomaly with the attraction of the terrain removed, and the subsequent restoration of these masses on the geoid as a thin layer according to the Helmert condensation method (Vanicek et al ., 1999), or some other method of restoration inside the geoid. We consider only the Helmert reduction, as this generally results in a small indirect effect, δN _TE The gravitational attraction of the topographic masses above the geoid may be derived from the potential given by Newton’s density integral and in spherical approximation is (Jekeli and Serpas, 2003) where the nominal value of the presumably constant mass density is ρ = 2670 kg/m ³ , and where r is the radial coordinate of the evaluation point, H is topographic height above the geoid, and If the terrain is radially condensed onto the geoid (sphere), thus creating a surface layer with density, ρH ( ξ ' ) , the corresponding gravitational attraction is given by The integrals in equations (5) and (8) are well defined if the evaluation points are not on the surface, r ≠ R + H ( ξ ) and r ≠ R, respectively, since then ℓ_h 0 = and ℓ ₀ = 0, if ψ = 0 . For topo δg _topo , the computation on the surface may be done numerically by omitting the integration near the evaluation point and substituting the formula for the attraction of a vertical prism of height, H . The layer effect is computed only for points with r > R , thus avoiding the singularity. The mathematical removal of the gravitational effect of the topography from the free-air anomaly results in the refined Bouguer anomaly, and the subsequent introduction of the equivalent mass layer yields the Helmert anomaly: There are two schools of thought on the requirement to apply a downward continuation of these anomalies to the geoid, either before or after adding the appropriate layer effect. That the Bouguer anomalies tend to be somewhat smoother argues in favor of their downward continuation, rather than the Helmert anomalies, which is also proposed here specifically for the airborne data. The process of combining airborne and terrestrial gravimetric data requires some type of downward continuation of the airborne data to the geoid. It is shown in Section 4 with simple numerical analyses that this combination is achieved favorably if both the airborne and the terrestrial data are first reduced to refined Bouguer anomalies. Indeed, the gravitational effect of the terrain is fundamentally different for both data types, where terrestrial data are affected by masses both below and above the measurement level, while airborne data are affected only by masses below. This difference creates a large bias between the corresponding free-air anomalies that is decreased significantly in the corresponding Bouguer anomalies. However, this also creates a complication since the downward continuation of Bouguer anomalies is potentially more difficult than free-air anomalies. For the terrestrial data, the downward continuation of Bouguer anomalies to the geoid is sometimes neglected as being a small effect (Jekeli and Serpas, 2003); however, see also (Huang and Véronneau, 2005). Certainly, it cannot be neglected for airborne data. One should also realize that downward continuation generally is possible by analytic methods only for harmonic fields, that is, in free space between level surfaces where Laplace’s field equation holds. Therefore, in any case, masses above the geoid should be removed prior to downward continuation to the geoid. Downward continuation of data is an unstable process since it is founded on the natural attenuation of the gravitational field in the opposite, upward direction. Hence, errors in the data at the high measurement surface are amplified when continued to the lower surface. The numerical procedure thus depends very strongly on proper filtering of the spectral content of the data, suppressing noise, but permitting sufficient signal content to be recovered at the lower surface (Novak and Heck, 2002). Several methods exist to downward-continue free-air gravity anomalies; they are based typically on inverting the field solution of a boundary-value problem. For example, Poisson’s integral (Hofmann-Wellenhof and Moritz, 2005) solves the Dirichlet problem and may also be formulated for gravity anomalies. Its inversion has the form where a r is the radial coordinate at altitude, and Δg ⁽⁰⁾ ( R, ξ ) = Δg ( r_a , ξ ) . The formulation is valid if no masses intervene between the surfaces, r = R and r = r_a . Another method is least-squares collocation (LSC, Moritz, 1980), which is also developed specifically for harmonic functions in free-space, and thus may be employed for harmonic continuation of functionals of the gravitational field. Another popular method is based on transforming equation (11) into the spectral domain, where the use of the Fast Fourier Transform (FFT) offers efficiency in computation (Bláha et al ., 1996; Hwang et al ., 2007). Our numerical experiments on simulated fields confirm that LSC and the Poisson inversion methods yield similar results under ideal conditions and are better than the FFT-based method. On the other hand, these same experiments also reveal that these methods applied to Bouguer anomalies fail significantly. More specifically, they fail to give accurate results for the downward continuation of the gravitational attraction of the terrain. It is conjectured that the reason for the failure is rooted in the fact that this attraction is not evaluated as a solution to a boundary-value problem; rather it is a forward computation based on Newton’s density integral. Of course, both ways of constructing a gravitational field are legitimate, either by means of a solution to a boundaryvalue problem, or directly modeling the masses that generate the field; and they should give the same field if there are no modeling errors. But the inversion methods and the forward modeling methods both have inherent approximations, and their numerical implementations make further different types of approximations. These practical inconsistencies conspire to make problematic the downward continuation of the forward-modeled δg _topo using a method based on inverting a solution to a boundary-value problem. The downward continuation of δg _topo should, therefore, also derive from a forward model. In fact, the radial derivatives, ∂^kδg _topo / ∂r^k , of the gravitational attraction at altitude due to the topographic masses may be computed similarly using forward models. Thus, a simple Taylor-series expansion to sufficiently high degree performs the required analytic continuation: It is noted that the downward continuation to the geoid is an analytic (harmonic) continuation that does not yield the actual gravitational attraction of the terrain at the geoid. However, this is of no consequence since the masses anyway have been removed for the Bouguer anomaly. No such straightforward continuation is possible for the free-air anomaly since it is modeled not with a forward model, but rather on the basis of a boundary-value problem. Therefore, to downward-continue the Bouguer anomaly, we may apply standard techniques (e.g., Poisson’s inverted integral) to the free-air anomaly component and the Taylor expansion to the terrain effect: Although separating the downward continuation in this way could pose theoretical problems if the Bouguer anomaly were a measured quantity, the procedure would be legitimate in this case since the Bouguer anomaly is constructed from two distinctly defined fields. Nevertheless, the downward continuation of δg _topo was not implemented for the present study. For South Korea, there is a considerable consolidated data base of 18677 terrestrial gravity anomalies (as of 2013), densely distributed in many parts of the country, but also sparsely distributed over other, mostly mountainous parts ( Figure 1 , left). Many of these data were collected since the 1990s by various institutions, including the Korea Institute of Geoscience and Mineral Resources (KIGAM), the Pusan National University (PNU), and the National Geographic Information Institute (NGII), among others. There is also a data set from a comprehensive airborne survey of gravity conducted over most of the country and its coastal regions in December 2008 and January 2009. These data were obtained along straight north-south tracks spaced at about 10 km, plus several east-west tracks spaced at about 50 km ( Figure 1 , right). The along-track spacing of the data is about 800 m, which, however, is higher spatial resolution than warranted due to an applied 150-s Butterworth filter. A representative along-track power spectral density (PSD) of the airborne data ( Figure 3 ) clearly shows the filter cutoff is near the frequency, f_c = 1×10 ⁻⁴ [cy/m] , or 10 km wavelength, corresponding to a resolution of 5 km. The terrestrial data, on the other hand, contain spectral information at much shorter wavelengths and with greater power at most frequencies, but their spatial resolution varies from tens of meters to tens of kilometers. Defining an equiangular 30"×30" geographic grid for South Korea, Table 1 quantifies the coverage of data, where a 30" cell is counted as contributing to the field representation if it contains at least one data point. At this resolution (about 1 km), the available terrestrial and airborne data cover only about 18% of the total land area of South Korea; and, the airborne data do not even have gravimetric information at this resolution. The situation is much improved for 2.5' resolution where these data cover about 84%. Nearly the entire country (97%) is represented at 5' resolution. Previous studies (Jekeli et al ., 2009) demonstrated that data resolution is a key factor in achieving high accuracy in the geoid height determination. It was predicted (ibid.) that for South Korea a 5' data resolution would achieve geoid height precision not better than 5-6 cm (standard deviation). It is a lower bound since it does not include errors coming from other sources, such as measurement error or errors in data reduction. Thus, although Table 1 indicates that 5' resolution has been achieved essentially with the combined terrestrial and airborne data, the accuracy of the value represented by the data in a 5' cell may vary considerably, especially if only a few data exist in that cell. With 2.5' resolution, the best achievable geoid precision was estimated (ibid.) to be about 2-2.5 cm (standard deviation). From Table 1 , it follows that some improvement beyond the 5 cm level may be possible. However, considering solely the airborne data, which have this resolution (along the track), we see that they cover only about 44% of the land area. Thus, it is imperative, firstly, that the terrestrial data are combined with the airborne data; and, secondly, that this combination is achieved with as much consistency as possible in order to avoid significant biases in the combined data set. Comparing the averaged terrestrial and airborne free-air anomalies in the common cells at any particular resolution, the statistics of their differences show a definite bias of about 8 mgal, as well a substantial standard deviation of 10-12 mgal ( Table 2 ). Understanding the source of this bias and correcting for it, therefore, is a major component in the data combination. Closer inspection (not shown) reveals that it is not a constant bias for all data, which is also seen indirectly in the PSDs ( Figure 2 ), where the terrestrial data generally are rougher in spectral content and have higher amplitudes at different spatial frequencies than the airborne data. If the gravitational attractions due to the topographic masses above the geoid are removed from the terrestrial and airborne data then the resulting refined Bouguer anomalies agree much better as shown in Table 3 . For each resolution, the Bouguer anomalies were computed at the given terrestrial and airborne data points using a 30" -resolution topographic model derived from the Shuttle Radar Topography Mission (SRTM, Farr et al ., 2007) and subsequently averaged over each of the common grid cells (as enumerated in Table 1 ). We see a practical elimination of the bias between the data sets and a significant reduction in the standard deviation. This then supports the proposition to combine the terrestrial and airborne data in terms of the refined Bouguer anomaly, as suggested in Section 3, rather than the free-air anomaly. The ultimate test of the quality of the gravity anomaly data is in the final desired product, the geoid height. The following setup was used to test the preliminary combination of airborne and terrestrial data. The data offset, N ₀ , in equation (2) was ignored as it cannot be determined from gravimetric data alone; however, it can be estimated in combination with GPS and leveling data, as shown below. The integration in this equation was limited to an area defined in latitude and longitude by 34° ≤ ϕ ≤ 38.5° and 125° ≤ λ ≤130° , respectively. This limitation in integration area is justified by the use of a high-degree reference model, given by EGM08 (Pavlis et al. , 2012) up to degree, n _max = 360 (30' resolution). Stokes’s function, S ( ψ ), was modified simply by removing its harmonic components with degrees less than or equal to M =120 according to the ideas developed by Wong and Gore (1969). It was found by simulation analyses that this is the optimal modification of this type for the given maximum resolution of the reference field. The full EGM08 model ( n _max = 2160) was used to compute free-air anomalies at all empty cells of an assumed data resolution, 30", 1' , or 2' . These include all otherwise empty cells on land, as indicated in Table 1 , and all ocean areas in the integration area. The Bouguer anomalies were then computed for all cells using the 30" SRTM terrain model. In cells that are common to the terrestrial and airborne data sets, the values were combined by simple averaging. Downward continuation was performed only on the airborne free-air anomaly, as an option, since the downward continuation of the Bouguer anomalies requires further study. Helmert anomalies were obtained finally using the same SRTM terrain model, according to equation (10). The indirect effect was computed using formulas given by Wichiencharoen (1982). We computed geoid heights with these data and procedures using equation (2) on a grid with the same resolution as the gravity anomaly. Their accuracy may be determined by comparing them to independently derived geoid heights, in this case, N obtained via equation (1), where the ellipsoidal height, h , was determined by GPS, and the orthometric height, H , was determined by precise leveling methods. It is estimated that the accuracy of these GPS/leveled geoid heights is about 2 cm (standard deviation). However, an additional small error enters when the gravimetric geoid heights are interpolated from the grid to the GPS/leveled point. Their point distribution, shown in Figure 3 , is part of the Unified Control Point (UCP) geodetic network newly established for South Korea with 10 km resolution. Table 4 summarizes the statistics of the differences at these points for different types and resolution of gravimetric data. As the data resolution decreases and the airborne data contribute more to the total field representation, the geoid precision based only on airborne data (with EGM08 fill-in) improves steadily. The opposite is evident for the terrestrial data – the geoid height precision deteriorates as less of the EGM08 contributes to the total data set. Combining the airborne and terrestrial data generally does not improve the geoid height values based only on airborne data, but yields more precise values than for the terrestrial-only case. Downward continuation of the airborne free-air anomaly offers a significant advantage, and applying the terrain effect also reduces the standard deviation in geoid height errors. It is anticipated that applying the downward continuation to the airborne gravitational attraction of the topographic masses could further reduce the geoid height errors. We note that the overall result for the geoid height obtained by Bae et al . (2012) is consistent with results reported in this paper. Their precision of 5.5 cm was achieved with a standard downward continuation of the airborne data by LSC. The mean errors in the geoid height are substantial and fairly consistent for different data resolutions if the terrain effect is included. This bias, in fact, is an estimate of the constant, N ₀ ,= i−n 1e4q.5u±at1io.8n c(m2). The values listed in Table 4 are consistent with N ₀ = −14.5 ±1.8 cm determined by Jekeli et al . (2012) using EGM08 and an older set of GPS/ leveled geoid heights. Further, it is evident in Table 4 that without the terrain reductions, and by implication without proper accounting of the biases between the terrestrial and airborne data, that the gravimetric geoid height may have a bias error of up to 10 cm. The determination of precise geoid heights requires high-resolution gravity data. The combined terrestrial and airborne data in South Korea equally contribute to about 2.5' (~5 km) resolution over 84% of the country and 5' (~10 km) resolution over the remaining parts. Preliminary analysis shows that geoid heights can be determined with precision better than 6 cm (standard deviation) when compared to more precise GPS/leveled geoid heights. An essential part in the combination of the terrestrial and airborne data is the gravitational attraction of the topographic masses above the geoid, which cause a significant, ~8 mgal, bias between the two data types. Accounting for this through the Bouguer reduction practically eliminates the bias and improves the geoid determination. Further improvements are anticipated and being investigated with proposed methods to apply more rigorous downward continuation to the geoid of the terrestrial and especially the airborne Bouguer anomalies.