O R I G I N A L R E S E A R C H P A P E R
Bathymetry enhancement by altimetry-derived gravity anomalies
in the East Sea (Sea of Japan)
Kwang Bae Kim•Yu-Shen Hsiao•Jeong Woo Kim• Bang Yong Lee•Yi Kyun Kwon•Chang Hwan Kim
Received: 10 February 2010 / Accepted: 18 November 2010 / Published online: 4 December 2010 Ó Springer Science+Business Media B.V. 2010
Abstract The gravity-geologic method (GGM) was used to enhance the bathymetry of the East Sea (Sea of Japan) with satellite altimetry-derived free-air gravity anomalies and shipborne depth measurements. By comparison with the bathymetry model of Smith and Sandwell’s (SAS) approach (1994), GGM was found to have an advantage with short wavelength (B12 km) components, while SAS better predicts longer wavelength (C25 km) components, despite its dependency on density contrast. To mitigate this limitation, a tuning density contrast of 10.25 g/cm3 between seawater and the seafloor was primarily estimated by the downward continuation method and then validated by the check points method with GGM. Similarly, SAS is limited by the ‘‘A’’ value in low-pass part of the Wiener filter, which defines the effective range of the wavelength components on bathymetry. As a final result, we present an enhanced GGM bathymetry model by integrating all available data.
Keywords Bathymetry enhancement Altimetry-derived free-air gravity anomalies East Sea (Sea of Japan) Gravity-geologic method (GGM) Density contrast
Introduction
Seafloor topography, or bathymetry, has traditionally been mapped by shipborne echo sounding measurements that are time-consuming and, because of the limited and uneven distribution of the measurements, often biased. However, the advent of the satellite radar altimeter has made it possible to estimate global bathymetry more economically and accurately. Since the late 1970 s, sev-eral satellite radar altimeters, such as Seasat (Born et al.
1979), Geosat (Cheney et al. 1986), ERS-1 and 2 (Got-tschalk 1991; Francis et al. 1995), and TOPEX/Poseidon (Fu et al. 1994) have provided dense sea surface topo-graphy by measuring the distance between the satellite and the sea surface.
These measurements can be converted to geoid, free-air gravity anomalies (FAGA) and bathymetry. Since the variations of FAGA are theoretically correlated with undulations of crustal density variations of the local bed-rock beneath the ocean floor, dense FAGA data that have been derived from satellite radar altimeter measurements can be used to predict bathymetry with the assumption that there is no horizontal density variation. Because of this limitation, shipborne bathymetry is required as a constraint. Several methodologies for recovering bathymetry have been studied in recent decades. Dixon et al. (1983) com-pared the predicted bathymetry by one-dimensional filter-ing of the geoid heights obtained from tracks of the Seasat satellite altimeter (Born et al. 1979) with observed K. B. Kim Y.-S. Hsiao J. W. Kim (&)
Department of Geomatics Engineering, University of Calgary, 2500 University Drive NW, Calgary, AB T2N 1N4, Canada e-mail: jw.kim@ucalgary.ca
Y.-S. Hsiao
Department of Civil Engineering, National Chiao Tung University, Hsinchu, Taiwan, ROC
B. Y. Lee
Korea Polar Research Institute, Incheon 406-840, Korea Y. K. Kwon
Korea Institute of Geoscience and Mineral Resources, Daejeon, 305-350, Korea
C. H. Kim
East Sea Research Institute, Korea Ocean Research and Development Institute, Uljin 767-813, Korea DOI 10.1007/s11001-010-9110-0
conventional bathymetry in the along-track directions. This study was followed by other one-dimensional along-track methods (Baudry et al. 1987; Craig and Sandwell 1988; Vogt and Jung 1991; Jung and Vogt 1992) and by two-dimensional bathymetry approaches (Baudry and Calmant
1991; Smith and Sandwell 1994; Calmant 1994; Sichoix and Bonneville 1996; Ramillien and Cazenave 1997; Hwang1999; Calmant et al.2002).
Calmant and Baudry (1996) described a comprehensive summary for bathymetry modeling of deep seafloor fea-tures, such as seamounts and ridges, on a global scale from satellite altimetry-derived FAGA. Bathymetric estimates derived from satellite altimetry and shipboard echo sound-ing data were reviewed by Sandwell and Smith (2001).
One of the useful techniques for estimating bathymetry is the gravity-geologic method (GGM). Originally, GGM was applied for mapping bedrock elevations beneath gla-cial drift deposits (Ibrahim and Hinze 1972; Adams and Hinze 1990; Nagarajan 1994). In terrestrial applications, GGM has been investigated for estimating the depth to bedrock from observed Bouguer gravity anomalies (BGA) and borehole measurements between the surface and the basement at the control points. One limitation in terrestrial application is that GGM predictions could be distorted because density varies between the glacial sediments and the bedrock. GGM has also been applied in marine settings, where the density contrast between seawater and the sea-floor bedrock is very homogeneous. The method can be used for predicting the ocean bottom topography by using FAGA and depth measurements instead of BGA and borehole measurements.
In the past 40 years, several research studies have been performed for bathymetry predictions using GGM from satellite altimetry-derived FAGA. GGM was adapted for bathymetric estimations in the Barents Sea and around Greenland (Roman1999). To implement GGM, a density contrast between seawater and the ocean bottom rock should be estimated. A density contrast of 9.0 g/cm3was estimated for GGM and used in the Drake Passage, Antarctica (Kim et al.2010). Determination of the effective density contrast estimated by the downward continuation method was recently implemented to predict bathymetry by GGM with an accuracy of less than 40 m (Hsiao et al.2010).
In this study, combined shipborne and satellite altimetry measurements were used to estimate a 1 9 1 arc-minute bathymetry model, by using dense altimetry-derived FAGA where sparse shipborne depth measurements were available. The study area (36–39°N and 130–133°E) is located in the East Sea (Sea of Japan) between Korea and Japan (Fig.1). As shown in Fig.1, the Ulleung Basin is located south of Ulleung and Dok Islands.
Data used in this study included shipborne depths and gravity data from the National Geophysical Data Center
(NGDC, http://www.ngdc.noaa.gov) GEODAS database (National Geophysical Data Center 2008), the National Oceanic and Atmospheric Administration (NOAA, http:// www.noaa.gov) and 1 9 1 arc-minute altimetry-derived FAGA (Sandwell and Smith 2009) generated by Geosat and ERS-1 satellite altimeters. Fig.2a shows the distribu-tion of the shipborne depths obtained from NGDC and the 1 9 1 arc-minute Earth topographical database 1 (ETO-PO1) model (Amante and Eakins2009) as the background. In the left panel of Fig.2a, white dots (18,317 points) have only shipborne depth measurements, and black dots (6,736 points) contain both depth and gravity measurements.
Figure2b shows the altimetry-derived FAGA in the study area, which was obtained from Scripps Institution of Oceanography (SIO, http://www.sio.ucsd.edu), University of California at San Diego (Sandwell and Smith2009). The original data were re-gridded to 1 9 1 arc-minutes by interpolation with the Generic Mapping Tools (GMT) software package (Wessel and Smith1998), because there were the distortions in the high latitudes that were made when the map was generated from the Mercator projection. To determine the density contrast, we first estimated a density contrast value by the downward continuation method and then evaluated to determine if the value was within the acceptable tuning density contrast range by GGM using check points (Hsiao et al. 2010). We also adopted Smith and Sandwell’s (SAS) method (1994) to estimate bathymetry from FAGA in this study. The two bathymetry models were compared with other bathymetry models, including the NGDC shipborne data and the ETOPO1 model. Furthermore, we validated our bathymetry models Fig. 1 Location map. Study area (36–39°N and 130–133°E) is delineated by a red box. The blue box is the boundary of the test area used for the density contrast determination. The 1 9 1 arc-minute bathymetry by ETOPO1 is superimposed as a background
with independently measured shipborne data by the Korea Institute of Geoscience and Mineral Resources (KIGAM,
http://www.kigam.re.kr) in April 2000. In this paper, we also present an enhanced bathymetry model by integrating all shipborne and satellite measurements in the study area as the final result of the study.
Bathymetry prediction by the gravity-geologic method Concept of the gravity-geologic method
GGM was originally developed to predict the depth-to-basement using a density contrast between glacial sedi-ments and bedrock from known gravity anomalies and
borehole measurements (Ibrahim and Hinze1972; Adams and Hinze 1990; Nagarajan1994).
In general, the observed BGA includes two gravity effects that are composed of the residual gravity generated from local bedrock variations and the regional gravity generated from deeper mass variations. In GGM applica-tions, we need the residual gravity effect, which yields short wavelength signals, to predict the final depth-to-bedrock. The regional gravity is computed by subtracting the residual gravity from the observed BGA if the depth measurements are known at the control points. However, at sites where the depth is not measured, we can estimate the residual gravity for predicting the bedrock undulations by subtracting the regional gravity, which is generated by gridding the regional gravity yielded at sites of measured Fig. 2 aThe 1 9 1 arc-minute
bathymetry by ETOPO1 with superimposed ship track measurements by the National Geophysical Data Center (NGDC). Black dots denote points with depth and gravity, while white dots have depth only. The test area for the determination of density contrast is also shown (with control points as white circles and check points as black triangles). b Satellite altimetry-derived free-air gravity anomalies (FAGA) (Sandwell and Smith2009) in the study area. Attributes listed for this and subsequent maps include the amplitude range (AR = minimum and maximum values), amplitude mean (AM), amplitude standard deviation (ASD), and amplitude unit (AU)
depth, j, from the observed BGA. The control points with measured depth are denoted with subscript, j, while the control points with unmeasured depth are indicated with subscript i in this study.
Local bedrock variations under sediments yield a shorter wavelength gravity field, while deeper mass variations generate a longer wavelength gravity field. These two gravity fields comprise the observed Bouguer gravity (gOBS) by combining the residual gravity (gRES) and the regional gravity (gREG) as in the following equation: gOBSðiÞ ¼ gRESðiÞ þ gREGðiÞ: ð1Þ
We assume that the measured depths at control points jk in Fig.3were used to estimate the residual bedrock gravity field (gRES) that generates the shorter wavelength effect from a simple Bouguer slab formula:
gRESðjÞ ¼ 2pGðDqÞðEðjÞ DÞ; ð2Þ
where G is the gravitational constant, 6.672 9 10-8cm3/ gsec2; Dq is the density contrast (g/cm3); E(j) is the bed-rock elevation at the jthcontrol point; and, D is the deepest depth of the control points as a reference datum. Both E(j) and D are measured in meters.
The residual gravity represents the effect of the bedrock surface at the control points and is removed from the observed Bouguer gravity to estimate the regional gravity field (gREG) that represents the longer wavelength effect: gREGðjÞ ¼ gOBSðjÞ gRESðjÞ: ð3Þ
The regional gravity, gREG(i), at sites of unmeasured depth, i, can be estimated by interpolation as a grid by using the generated regional gravity, gREG(j), at sites of measured depth, j. The estimated regional gravity, gREG(i), at site i is removed from the observed gravity, gOBS(i), to estimate the residual gravity, gRES(i):
gRESðiÞ ¼ gOBSðiÞ gREGðiÞ: ð4Þ
By rearranging Eq.2, the elevation of the bedrock above the datum where bedrock depths are unmeasured at site i can be estimated with the following formula:
EðiÞ ¼ gRESðiÞ
2pGDqþ D: ð5Þ
The estimation of bedrock topography using GGM is influenced by the distribution and number of the depths and gravity at the control points, the observed gravity, and the density contrast between the bedrock and the overlying sediments. Where the shipborne depth measurements are sparse over oceanic areas, GGM can be used to predict bathymetry from satellite altimetry-derived dense FAGA. GGM applications in ocean areas use shipborne depths and FAGA, instead of bedrock elevation and Bouguer gravity at control points, in order to predict the ocean bottom topography under seawater (Roman1999). FAGA data are approximately at the geoid in the ocean and, thus, do not need correction for height variation and Bouguer mass effects.
The density contrast between seawater and the seafloor varies as a result of changes in mass under the ocean bottom. In GGM applications, a geologically reasonable density contrast between seawater (1.03 g/cm3) and the seafloor bedrock (2.70 g/cm3) is assumed to be 1.67 g/cm3. The density of the seafloor bedrock in the crust is assumed to be the mean value between 2.67 g/cm3 and 2.73 g/cm3 (Jin
1995; Kim et al. 2010). In reality, the density contrast is much larger than the theoretical value (1.67 g/cm3), because satellite altimetry-derived FAGA compiled on the sea sur-face in GGM are caused by the sources at the ocean bottom (Roman1999; Kim et al.2010). For more effective bathy-metric prediction, a tuning value of density contrast that was determined from the control points of measured shipborne depths was used (Kim et al.2010; Hsiao et al.2010). Determination of density contrast
Check points with the gravity-geologic method
Because the shipborne depth measurements shown in Fig.2a have large gaps, we selected the southeast corner (36.0–37.46°N and 132.14–133.0°E) as a test area for estimating the density contrast. This test area is delineated by the blue box in Figs.1and2a, where the distribution of shipborne measurements is denser and more even. The right-hand panel in Fig. 2a shows control and check points within the test area with the ETOPO1 bathymetry model as the background. Of all 1,493 shipborne depth measure-ments within the test area, two thirds (992) of the points were used as control points (white circles). The remaining third (501) were used as check points (black triangles). Every third point of shipborne depths along shiptracks was Fig. 3 Geometry of the gravity-geologic method and the downward
selected as a check point to evenly distribute the control points within the test area. The control points in the test area were used to effectively determine the density contrast of the study area.
Figure4 is a trade-off diagram, which shows a root-mean-square (RMS) difference, its rate of change and a correlation coefficient between the control and check points. This trade-off diagram was developed by the check points method for picking an effective density contrast for GGM bathymetry predictions. If the selected density con-trast is too small, the GGM-predicted bathymetry would be very large. In contrast, if the selected density contrast is too big, the GGM-predicted bathymetry would be too close to ‘‘level n’’ (Fig.3). A range of density contrasts that sta-bilize the RMS difference by using the check points method was selected in the trade-off diagram. From the three curves in the trade-off diagram, we selected the range of acceptable density contrasts as 9.0 g/cm3 and higher, which is the shaded area in Fig.4.
The determination of density contrast between seawater and the ocean bottom topographic mass is an important factor for GGM bathymetry. We know that the density contrast (Dq) in Eq.5 controls the amplitude of bathy-metric estimates. Although it is not geologically reason-able, the check points were used to select an optimal density contrast in the GGM application (Nagarajan1994; Kim et al.2010).
Downward continuation method
In this study, the downward continuation method was used to select an effective density contrast within the acceptable range, i.e. 9.0 g/cm3and higher, determined by the check points method described in the previous section. Upward continuation is regarded as a smooth operation (Tikhonov
and Arsenin 1977), which is a well-posed problem whose solution will be unique; whereas, downward continuation is characterized as an unstable solution, which is an ill-posed problem whose solution will be non-unique and unstable. Upward continuation in the frequency domain from the gravity field at z = h1 plane to z = h2 plane can be explained as (Parker1977):
Gh2ðu; vÞ ¼ expð2pkdÞGh1ðu; vÞ; ð6Þ
where Gh1ðu; vÞ and Gh2ðu; vÞ are the two-dimensional
Fourier transforms of the gravity field at h1 and h2, respectively; u and v are the frequencies for x and y directions, respectively; k¼pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiu2þ v2; and d = h2- h1.
Downward continuation of the gravity field in the Fourier transform from z = h2plane to z = h1plane, by contrast, is given by inverting Eq.6 as:
Gh1ðu; vÞ ¼ expð2pkdÞGh2ðu; vÞ: ð7Þ
The sea surface gravity field using the downward continuation method was calculated for several downward continued levels, as shown in Fig.3, by the Gaussian filter in Eq. 7:
Gðu; vÞ ¼ Fðu; vÞGðu; vÞ expð2pkdÞ; ð8Þ where G represents the downward continued gravity field, G represents the original gravity field, d denotes the dis-tance of the downward continuation, and F is a Gaussian filter.
The effective density contrast for considering mass variations in the deepest ocean bottom was determined by the downward continuation method until reaching the seafloor topography. The density contrast from the esti-mated gravity ratio between the sea surface level and each downward continued level was computed, as shown in Fig.3. The density contrast (DqDWC) at point i of every
Fig. 4 Trade-off diagram for selecting the tuning density contrast range of 9.0 g/cm3and higher by the gravity-geologic method and a final density contrast of 10.25 g/cm3by the downward continuation method
downward continued level can be estimated as a multi-plying gravity ratio, DgDWC(i)/DgSEA(i), at point i of each downward continued level and sea surface level by DqSEA, which is the density contrast at the sea surface level as a constant of 1.03 g/cm3(Strykowski et al.2005). The den-sity contrast in each downward continued level was com-puted as the mean value of DqDWC.
Figure5shows the predicted density contrasts for every level at intervals of 200 m below the sea surface, along with the gravity ratios between the sea surface and the downward continued level. The density contrast of each level was calculated until the levels reached 2,037 m, which is the greatest depth in the test area. As shown in Fig.5, both the gravity ratio and the predicted density contrast gradually increased with the ocean depth.
From the two curves, 10.25 g/cm3, which is the pre-dicted density contrast at the deepest depth, was selected as a tuning density contrast that would be applied to the GGM bathymetry predictions, because the tuning density contrast is within the acceptable range determined by the check points method with GGM.
Bathymetry prediction by Smith and Sandwell’s approach
Review of Smith and Sandwell’s approach
The Parker formula (Parker 1972) can be described as a forward model used to convert seafloor topography (h) to gravity anomaly (Dg) in the two-dimensional Fourier domain: =½Dg ¼ 2pGðDqÞ expð2pkdÞX 1 n¼1 ð2pkÞn1 n! =½h n; ð9Þ
where G is the gravitational constant; Dq is the density contrast between seawater and the seafloor topography; d is the mean depth beneath the ocean surface; k is the radial frequency ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffik2
xþ ky2
q
; k is (kx, ky) = (1/kx, 1/ky), where (kx, ky) and (kx,ky) are frequencies and wavelengths in the x and y directions, respectively; and,=½is the Fourier transform operator.
If the magnitude of the seafloor topography undulation is much less than the mean ocean depth (d), the first term dominates because the series of non-linear terms of Eq.9
converge fast. The first term in Eq.9is only dependent on the mean depth below the ocean surface and is a linear relation-ship between the gravity anomaly and the seafloor topogra-phy in the Fourier domain as in the following equation: GðkÞ ¼ ZðkÞHðkÞ ¼ 2pGðDqÞ expð2pkdÞHðkÞ; ð10Þ where G(k) and H(k) are the Fourier transforms of the gravity anomalies (Dg) and the seafloor topography (h), respectively; Z(k) is the admittance function as a transfer function; and, k is the frequency in two-dimensions. Z(k) is a transfer function that is isotropic and spatially invariant. In the Z(k) terms of Eq.10, the 2pGDq term is a Bouguer constant as an infinite slab of material; and, the exp(-2pkd) term represents upward continuation (from seafloor to sea surface) through the mean ocean depth as exponential decay with increasing frequency.
Next, we estimated the seafloor topography from satel-lite altimetry-derived gravity anomalies using the inverse model of Eq.10. The inversion model is shown as the following equation in the Fourier domain:
HðkÞ ¼ Z1ðkÞGðkÞ ¼ 1
2pGDq expð2pkdÞGðkÞ: ð11Þ Smith and Sandwell (1994) (SAS) suggested that, in a limited band of intermediate wavelengths between about 15 and 160 km, the satellite altimetry-derived FAGA data have a correlation between gravity and bathymetry, because they can be implemented by attenuation of gravity by upward continuation in shorter wavelengths and by isostatic compensation of the seafloor topography at longer wavelengths. For a stable condition of downward continuation, which is represented as the term exp(2pkd) in Eq. 11, to estimate bathymetry from gravity anomalies, the transfer function Z-1(k) is suppressed by windowing the limited range of wavelengths of topography. Thus, a band-pass filter was used. The purpose of this filter is to stabilize the downward continuation in predictions of the seafloor topography (Smith and Sandwell 1994; Calmant and Baudry1996; Watts et al.2006).
Fig. 5 Density contrast at every level and gravity ratio between the sea surface and the downward continued level in the test area
In the 1994 SAS approach, the altimetry-derived gravity anomalies in the frequency domain, G0(k), are first band-pass filtered and downward continued to a mean depth (d) by:
GðkÞ ¼ G0ðkÞWðkÞ expð2pkdÞ; ð12Þ
where G(k) is the band-pass filtered and downward continued gravity field, and W(k) is the Wiener filter. The gridded shipborne depth measurements in the frequency domain, B0(k), are also band-pass filtered by the Wiener filter:
HðkÞ ¼ B0ðkÞWðkÞ; ð13Þ
where H(k) is the seafloor topography field.
In Eqs.12and13, the Wiener filter is composed of two filters, W1(k) as the high-pass filter and W2(k) as the low-pass filter:
WðkÞ ¼ W1ðkÞ W2ðkÞ; ð14Þ
where W1ðkÞ ¼ 1 exp½2ðpksÞ2, W2ðkÞ ¼ ½1 þ
Ak4expð4pkdÞ1, and s is the Wiener filter parameter. The ‘‘A’’ value chosen in W2(k) affects the cut-off wavelength scale at the resolution of the predicted bathymetry. Thus, the choice of the ‘‘A’’ value must be considered rigorously and should represent an appropriate compromise between noise reduction and signal detail. Unlike SAS (1994), who used 9,500 km4for the ‘‘A’’ in the low-pass filtering (i.e. W2(k)), we selected 6,233 km4in our study with a new gravity model (personal communication with David Sandwell, 2009). The gridded shipborne depth measurements in the frequency domain, B0(k), are then low-pass filtered by:
DLPðkÞ ¼ B0ðkÞ½1 W1ðkÞ; ð15Þ
where DLP(k) is the low-pass filtered gridded shipborne depth in the frequency domain.
Finally, the SAS bathymetry is estimated by:
SAS¼ dLPðxÞ þ S gðxÞ; ð16Þ
where dLP(x) and g(x) are the spatial domains of DLP(k) in Eq.15 and G(k) in Eq. 12, respectively; S is a scaling factor determined by linear regression between G(k) in Eq.12and H(k) in Eq.13; and, dLP(x) and Sg(x) represent the longer and band-limited wavelength bathymetry, respectively.
Results from Smith and Sandwell’s approach
The SAS bathymetry model is shown in Fig.6b. In SAS, the linear relationship between gravity anomalies and the seafloor topography was fundamentally used in the Fourier domain. Band-pass filtering was performed in the down-ward continued satellite altimetry gravity anomalies and the gridded shipborne depths using Eqs.12 and 13,
respectively, with a mean water depth of 1,868 m in the study area. g(x) in the spatial domain obtained from Eq.12
was converted into the seafloor topography by multiplica-tion with the scaling factor, S, which was empirically computed from linear regressions of band-pass filtered and downward continued gravity anomalies, g(x), and band-pass filtered gridded shipborne depths, h(x), in the study area.
The scaling factor of 2.20 m/mGal was used to convert band-pass filtered and downward continued altimetry-derived gravity anomalies into bathymetry. For longer wavelengths with regional effect, low-pass filtering of the gridded shipborne depth measurements was implemented in Eq.15. Finally, the SAS bathymetry model was gener-ated from two bathymetries at different wavelengths.
Comparison of bathymetry models
The 1 9 1 arc-minute bathymetry models by GGM, SAS, NGDC, and ETOPO1 are shown in Fig. 6. The NGDC gridded model was produced by the harmonic surface technique in GMT software (Wessel and Smith 1998). Statistical comparisons are summarized in Table1. The GGM model showed correlation coefficients of 0.998, 0.999, and 0.957 with SAS, NGDC, and ETOPO1 models, respectively. SAS also had a maximum correlation coeffi-cient of 0.998 with NGDC.
The reason why both GGM and SAS were most strongly correlated with the NGDC gridded model is that NGDC shipborne depths were utilized to estimate both bathymetry models. According to the comparison in Table 1, the RMS of the GGM model is smaller than those of the SAS, NGDC and ETOPO1 models.
Maps of point-by-point differences between the models are represented in Fig. 7, and Table2 summarizes their statistics. In Fig.7b, the difference between the GGM and NGDC gridded models (GGM–NGDC) represents good agreement with the smallest RMS (34.1 m) in Table2. However, the differences in shallow water around Ulleung Island in Fig.7c, e, f reached approximately -2,000 m. These large differences were caused by the steep bathym-etry gradient around the island.
The power spectral density of the bathymetry models was analyzed in the radial direction to examine energy variations in different wavelengths. In Fig. 8, the amplitude of the power spectral density along the y-axis is repre-sented as a unit of dB that is computed from 10log10(P), where P is a power. The wavelength along the x-axis is considered to range up to 300 km due to the size of the study area. As illustrated in Fig.8, the SAS model showed a high power spectral density for wavelengths longer than
25 km and presented low energy in short wavelengths less than roughly 12 km.
Both the GGM and SAS models showed similar energy levels in the medium wavelengths between 12 km and 25 km; however, the GGM model had more power in very
short wavelengths of less than approximately 12 km. These results may indicate that the GGM has better bathymetry in the short wavelength components; whereas, the SAS can effectively yield bathymetry in the medium to longer wavelength components.
Fig. 6 Bathymetry models by athe gravity-geologic method (GGM), b the Smith and Sandwell (SAS) approach (1994), c the National Geophysical Data Center (NGDC) and d Earth topographical database 1 (ETOPO1). Bathymetry estimates along Profiles #1 and #2 are contrasted. Red dots denote shipborne measurements with anomalous depths
Table 1 Statistic comparisons between bathymetry models by the gravity-geologic method (GGM), Smith and Sandwell’s (SAS) approach (1994), National Geophysical Data Center (NGDC), and Earth topographical database 1 (ETOPO1) in the study area (unit: meter)
Min. Max. Mean SD RMS Correlation coefficient (CC)
GGM SAS NGDC ETOPO1
GGM -3998.1 -49.8 -1853.4 693.7 1979.0 1 0.998 0.999 0.957
SAS -3238.8 0.0 -1853.4 720.7 1988.6 0.998 1 0.998 0.959
NGDC -3238.6 -52.8 -1857.7 693.4 1982.9 0.999 0.998 1 0.950
To validate the result of the power spectral density analysis, we directly compared the GGM and SAS models with NGDC shipborne measurements possibly with short
wavelength components. Before the comparison, the two models were interpolated into the NGDC shipborne loca-tions, which are shown as black dots in Fig.2a. Figure9a, Fig. 7 The differences between
four bathymetry models: GGM, SAS, NGDC and ETOPO1. aGGM – SAS, b GGM – NGDC, c GGM – ETOPO1, dSAS – NGDC, e SAS – ETOPO1, and f NGDC – ETOPO1
b show the absolute values of the differences between the NGDC ship depths and the GGM bathymetry and between the NGDC ship depths and the SAS bathymetry on the NGDC shipborne tracks.
As shown in Fig.9a, the minimum, maximum, mean, and standard deviation of the absolute differences are 0.0, 1222.7, 16.1, and 37.8 m, respectively. Ninety-three per-cent of the absolute differences were less than 37.8 m, which is the standard deviation. Only 1.5% of the absolute differences were more than 100 m in the NGDC ship tracks in the study area. In Fig.9b, the statistics of the absolute differences between the SAS bathymetry and the NGDC bathymetry on the NGDC ship tracks were 0.0, 656.7, 46.9, and 44.4 m for the minimum, maximum, mean, and stan-dard deviation, respectively. The absolute difference below 44.4 m (the standard deviation) was 60.5% and 10.9% for more than 100 m.
From the histograms of two absolute differences from the NGDC ship tracks and the two competing models presented in Fig.9c, d, we know that 95% of the absolute differences between NGDC and GGM were below 50.0 m; whereas, 67.4% of the absolute differences between NGDC and SAS were less than 50.0 m. Based on the small dif-ferences between the GGM model and the NGDC ship-borne depths along the NGDC ship tracks, we concluded that GGM is more suitable than the SAS method for short wavelength bathymetry prediction.
GGM bathymetry estimated using a tuning density contrast (10.25 g/cm3) and determined by the downward continuation method is enhanced with an accuracy of 28.4 m in the standard deviation (693.7 m in Table1), when comparing the standard deviation (722.1 m) of the GGM bathymetry predicted by the geologically reasonable density contrast (1.67 g/cm3) between seawater and the seafloor topography.
In addition to the NGDC shipborne depths, shipborne data measured in April 2000 by the Korea Institute of Geoscience and Mineral Resources (KIGAM), which are not included in the NGDC, NOAA databases, were com-pared with the GGM, SAS, and NGDC models along two profiles: Profile #1 (east–west direction) and Profile #2 (north–south direction), shown in white lines in Fig.6a, b, c. The results are shown in Fig.10.
The KIGAM shipborne depths in Profile #1 were strongly correlated with GGM, SAS, and NGDC with correlation coefficients of 0.982, 0.976, and 0.979, respectively. For Profile #2, the correlation coefficients were also very strong with 0.997, 0.998, and 0.998 for GGM, SAS, and NGDC, respectively. Table3summarizes the statistics of KIGAM, GGM, SAS, and NGDC depths along Profiles #1 and #2 of the KIGAM ship tracks.
Table4presents the statistics of the differences between KIGAM and GGM, SAS, and NGDC bathymetries along Profiles #1 and #2. In Profile #1, the RMS of the differ-ences between the KIGAM ship bathymetry and each of the GGM, SAS, and NGDC bathymetries were 43.5, 79.7 and 46.8 m, respectively; and, in Profile #2, the difference between KIGAM and GGM had the smallest RMS value (28.7 m).
However, because of the overestimated SAS bathymetry in shallow water areas in Fig.10a, b, the profiles between KIGAM and SAS had a discrepancy of about 100 m between 130.45 and 130.65°E in Profile #1 and between 37.65 and 37.72°N in Profile #2. These differences were caused by the scaling factor (S) in the SAS bathymetry model; therefore, we can decrease this difference by adjusting S. We know that short wavelength bathymetry is well recovered by the GGM model, because the GGM bathymetry indicated small RMS differences compared with the two profiles (Profiles #1 and #2) collected along Table 2 Statistical comparison of the differences between four
bathymetry models: GGM, SAS, NGDC, and ETOPO1 (unit: meter)
Min. Max. Mean SD RMS
GGM – SAS -1665.3 304.6 0.0 52.8 52.8 GGM – NGDC -1596.2 289.1 4.3 33.8 34.1 GGM – ETOPO1 -1877.7 967.3 -4.0 206.6 206.6 SAS – NGDC -119.5 372.4 4.3 53.6 53.8 SAS – ETOPO1 -1801.1 934.6 -4.1 205.8 205.9 NGDC – ETOPO1 -2128.4 946.0 -8.4 221.3 221.5
Fig. 8 The comparison of power spectral densities between GGM (red circles) and SAS (blue circles) bathymetry models. The unit of power spectral density is dB, which is computed from 10log(P), where P is a power
the KIGAM ship tracks with very short wavelength bathymetry.
After examination of the maps in Fig.6, we found that the three bathymetry models (GGM, SAS, and NGDC) had ship tracks with outlier data, when compared with the ETOPO1 model. The outliers were considered as anoma-lous depths and are shown as red dots in Fig.6a–d. We eliminated ship tracks with outliers, as shown in red dots in the left panel of Fig.2a. The eliminated shipborne depths in Fig.2a were not used for predicting an enhanced bathymetry model.
As a final result, we present an enhanced bathymetry model by using all shipborne measurements, including
KIGAM data, in GGM, with a density contrast of 10.25 g/ cm3 estimated by GGM and the downward continuation method in the test area (Fig.11a). The ETOPO1 model was also applied to the islands, where the topography is above sea level. The enhanced bathymetry model had a minimum, maximum, mean, and standard deviation of -3182.3, 756.0, -1858.9, and 699.7 m, respectively.
The white color in the new model shown in Fig. 11a denotes depth at or above sea level. Figure11b shows the difference between the enhanced and original GGM mod-els. Overall, the enhanced model was similar to the original GGM model in Fig. 6a, except for the shallow water around Ulleung and Dok Islands, where the new model Fig. 9 The absolute differences of bathymetry on the NGDC shipborne tracks a between NGDC and GGM and b between NGDC and SAS. The histograms of the absolute differences on the NGDC shipborne tracks c NGDC and GGM and d NGDC and SAS
provided a better estimation by ETOPO1 for the area above sea level. The correlation coefficient between the enhanced and original GGM models was 0.966.
Results and conclusions
In this study, the gravity-geologic method (GGM) was used to enhance a bathymetry model from NGDC shipborne depths and satellite altimetry-derived FAGA (Sandwell and Smith2009) in the area of 36–39° N and 130–133° E in the
East Sea. The denser shipborne measurements of the test area in the southeast part of the study area were used to calculate the density contrasts between the seawater and the ocean bottom topography. A tuning density contrast of 10.25 g/cm3, which was estimated by the downward con-tinuation method, was found to be in the acceptable range (9.0 g/cm3 and higher), as determined with control and check points.
According to the flexural isostatic compensation theory, the relationship between bathymetry and downward con-tinued gravity anomalies is correlated in the 15–160 km wavelength band (Smith and Sandwell 1994). We imple-mented the Smith and Sandwell’s (SAS) approach (1994) by using the linear relationship between the seafloor topography and the downward continued gravity anomalies in the study area. The GGM bathymetry had correlation coefficients of 0.998, 0.999, and 0.957 with the SAS, NGDC, and ETOPO1 bathymetry models, respectively. The difference between the GGM and NGDC gridded models from the results of comparisons of the four bathymetry models (GGM, SAS, NGDC gridded, and ETOPO1) presented smaller RMS than the differences with the others.
In addition, the standard deviation of the absolute dif-ferences between the GGM bathymetry model and the NGDC shipborne depths on NGDC ship tracks was below 37.8 m. From the power spectral density analysis, we showed that GGM has an advantage with short wavelength (B12 km) components, while the SAS provides a better prediction of longer wavelength (C25 km) components. In the medium wavelength range, the two methods perform similarly.
An enhanced GGM bathymetry model was presented as a final result by combining all available data.
We conclude that GGM better predicts bathymetry over short wavelengths (B12 km) than SAS. Moreover, the bathymetry calculated by GGM is easily upgradable with Fig. 10 Bathymetry estimate comparisons for Profiles #1 (a) and #2 (b), respectively, as shown in Fig.6. CC is correlation coefficient
Table 3 Statistic comparisons between the Korea Institute of Geo-science and Mineral Resources (KIGAM), GGM, SAS, and NGDC bathymetry along Profiles #1 and #2 (unit: meter)
Min. Max. Mean SD RMS
KIGAM Profile #1 -2032.5 -1225.2 -1487.5 185.9 1499.0 Profile #2 -2080.8 -1235.4 -1611.5 323.2 1643.4 GGM Profile #1 -1908.7 -1222.9 -1461.8 185.2 1473.4 Profile #2 -2078.0 -1221.1 -1610.0 309.7 1639.3 SAS Profile #1 -1949.9 -1133.7 -1442.1 233.0 1460.7 Profile #2 -2087.9 -1145.9 -1563.1 341.4 1599.7 NGDC Profile #1 -1931.1 -1216.5 -1464.1 197.0 1477.2 Profile #2 -2070.3 -1214.8 -1597.5 306.6 1626.5
Table 4 Statistics of the differences between KIGAM, GGM, SAS, and NGDC bathymetry in Profiles #1 and #2 (unit: meter)
Min. Max. Mean SD RMS KIGAM – GGM Profile #1 -123.8 74.6 -25.6 35.2 43.5
Profile #2 -62.7 57.4 -1.5 28.8 28.7 KIGAM – SAS Profile #1 -137.4 103.5 -45.4 65.8 79.7 Profile #2 -94.4 7.1 -48.5 27.2 55.5 KIGAM – NGDC Profile #1 -101.4 94.1 -23.4 40.7 46.8 Profile #2 -73.5 38.5 -14.0 26.7 30.1
future satellite and shipborne measurements. Because the accuracy of GGM is dependant on the density contrast, we present a method for estimating a suitable tuning density contrast by combining the GGM and downward continua-tion methods. The use of the Bouguer slab approximacontinua-tion
mitigated by establishing tuning densities stabilized the GGM predictions (Kim et al. 2010), especially for mod-eling the gravity effects of variable density and rugged bathymetric relief at distances up to several kilometers.
In fact, the estimated tuning density of 10.25 g/cm3was found to have a more analytical than physical significance, because it essentially mitigates limitations in using a sin-gle-density Bouguer slab for analyzing multi-density ter-rain with rugged relief components that require gravity terrain corrections (Nagarajan1994; Kim et al.2010). The GGM bathymetry predicted by the tuning density contrast (10.25 g/cm3) presents an improvement of 28.4 m over the GGM bathymetry estimated by the reference density con-trast (1.67 g/cm3) between seawater and the ocean bottom topography, although the tuning density is larger than a geologically reasonable value. The GGM model using the geologically reasonable density contrast is not realistic. It implies that the influence of density contrast for GGM is unavoidable. Because the SAS approach includes a band-pass filter, it may result in underestimating the short wavelength components.
Geology and bathymetry of the East Sea recently became an important issue, not only because understanding the origin of the basins in the East Sea is important for evaluation of natural resources (i.e. gas hydrate and hydrocarbon), but also in terms of exclusive economic zones (EEZs) and territorial boundaries between the neighboring countries. As a result, the GGM method is a useful and convenient approach for updating bathymetry when new satellite altimetry mission gravity data and shipborne depths become available.
Acknowledgments This study was supported by the COMPAC Project (PE10030) and the Korea Arctic Multidisciplinary Program (KAMP, PP10090) of the Korea Polar Research Institute.
References
Adams JM, Hinze WJ (1990) The gravity-geologic technique for mapping varied bedrock topography. In: Ward SH (ed) Geo-technical and environmental geophysics. Environmental and groundwater, vol III. Society of Exploration Geophysicists, Tulsa, OK, pp 99–106
Amante C, Eakins BW (2009) ETOPO1 1 arc-minute global relief model: procedures, data sources and analysis. NOAA Technical Memorandum NESDID NGDC-24
Baudry N, Calmant S (1991) 3-D modelling of seamount topography from satellite altimetry. Geophys Res Lett 18(6):1143–1146 Baudry N, Diament M, Albouy Y (1987) Precise location of
unsurveyed seamounts in the Austral archipelago area using SEASAT data. Geophys J R Astron Soc 89:869–888
Born GH, Dunne JA, Lamb DB (1979) Reports: Seasat mission overview. Science 204:1405–1406
Calmant S (1994) Seamount topography by least-squares inversion of altimetric geoid heights and shipborne profiles of bathymetry and/or gravity anomalies. Geophys J Int 119:428–452
Fig. 11 a Enhanced bathymetry model and b differenced (Fig.11a - Fig. 6a) GGM bathymetry models. Note that the two color scales are different
Calmant S, Baudry N (1996) Modeling bathymetry by inverting satellite altimetry data: a review. Mar Geophy Res 18:123–134 Calmant S, Berge-Nguyen M, Cazenave A (2002) Global seafloor topography from a least-squares inversion of altimetry-based high-resolution mean sea surface and shipboard soundings. Geophys J Int 151:795–808
Cheney R, Douglas B, Green R, Miller L, Milbert D, Porter D (1986) The GEOSAT altimeter mission: a milestone in satellite oceanography. EOS Trans Am Geophys Un 67(48):1354–1355 Craig CH, Sandwell DT (1988) Global distribution of seamounts from
Seasat profiles. J Geophys Res 93(B9):10408–10420
Dixon TH, Naraghi M, McNutt MK, Smith SM (1983) Bathymetric prediction from SEASAT altimeter data. J Geophys Res 88(C3): 1563–1571
Francis CR, Graf G, Edwards PG, McCaig M, McCarthy C, Lefebvre A, Pieper B, Pouvreau PY, Wall R, Weschler F, Louet J, Schumann W, Zobl R (1995) The ERS-2 spacecraft and its payload. ESA Bull 83:13–31
Fu LL, Christensen EJ, Yamarone CA Jr, Lefebvre M, Menard Y, Dorrer M, Escudier P (1994) TOPEX/POSEIDON mission overview. J Geophys Res 99(C12):24369–24381
Gottschalk D (1991) ERS-1 mission and system overview. Die Geowissenschaften 9:100–101
Hsiao YS, Kim JW, Kim KB, Lee BY, Hwang C (2010) Bathymetry estimation by gravity-geologic method: investigation of density contrast predicted by downward continuation. Terr Atmos Ocean Sci (accepted for publication)
Hwang C (1999) A bathymetric model for the South China Sea from satellite altimetry and depth data. Mar Geod 22:37–51 Ibrahim A, Hinze WJ (1972) Mapping buried bedrock topography
with gravity. Ground Water 10(3):18–23
Jin YK (1995) Crustal structure of the South Shetland trench and the Shackleton fracture zone off the northern Antarctic Peninsula. Ph.D. dissertation, Seoul National University
Jung WY, Vogt ER (1992) Predicting bathymetry from Geosat-ERM and shipborne profiles in the South Atlantic Ocean. Tectono-physics 210:235–253
Kim JW, von Frese RRB, Lee BY, Roman DR, Doh SJ (2010) Altimetry-derived gravity predictions of bathymetry by gravity-geologic method. Pure Appl Geophys. doi: 10.1007/s00024-010-0170-5(accepted for publication)
Nagarajan R (1994) Gravity-geologic investigation of buried bedrock topography in northwestern Ohio. M.Sc. thesis, Department of Geological Sciences, Ohio State University
National Geophysical Data Center (2008) http://www.ngdc. noaa.gov/mgg/gdas/ims/trk_cri.html. US Department of Com-merce, National Oceanic and Atmospheric Administration, National Geophysical Data Center, Boulder, Colorado, USA Parker RL (1972) The rapid calculation of potential anomalies.
Geophys JR Astron Soc 31:447–455
Parker RL (1977) Understanding Inverse Theory. Ann Rev Earth Planet Sci 5:35–64
Ramillien G, Cazenave A (1997) Global bathymetry derived from altimeter data of the ERS-1 geodetic mission. J Geodyn 23(2):129–149
Roman DR (1999) An integrated geophysical investigation of Greenland’s tectonic history. Ph.D. dissertation, Department of Geological Sciences, Ohio State University
Sandwell DT, Smith WHF (2001) Bathymetric estimation. In: Fu LL, Cazenave A (eds) Satellite altimetry and earth sciences. A handbook of techniques and applications. Academic press, San Diego, pp 441–457
Sandwell DT, Smith WHF (2009) Global marine gravity from retracked Geosat and ERS-1 altimetry: Ridge segmentation versus spreading rate. J Geophys Res 114(B01411). doi: 10.1029/2008JB006008
Sichoix L, Bonneville A (1996) Prediction of bathymetry in French Polynesia constrained by shipboard data. Geophys Res Lett 23(18):2469–2472
Smith WHF, Sandwell DT (1994) Bathymetric prediction from dense satellite altimetry and sparse shipboard bathymetry. J Geophys Res 99(B11):21803–21824. doi:10.1029/94JB00988
Strykowski G, Boschetti F, Papp G (2005) Estimation of the mass density contrast and the 3D geometrical shape of the source bodies in the Yilgarn area, Eastern Goldfields, Western Austra-lia. J Geodyn 39(5):444–460
Tikhonov AN, Arsenin VY (1977) Solutions of Ill-posed problems. Winston-Wiley, New York
Vogt ER, Jung WY (1991) Satellite radar altimetry aids seafloor mapping. EOS Trans Am Geophys Un 72(465):468–469 Watts AB, Sandwell DT, Smith WHF, Wessel P (2006) Global
gravity, bathymetry, and the distribution of submarine volcanism through space and time. J Geophys Res 111:B08408. doi: 10.1029/2005JB004083
Wessel P, Smith WHF (1998) New improved version of the Generic Mapping Tools released. EOS Trans Am Geophys Un 79(579). doi:10.1029/98EO00426
