Gravity Anomalies from retracked ErS and Geosat Altimetry over the Great Lakes : Accuracy Assessment and Problems

Lake waters are a useful area for investigating the response of the Earth’s gravity field to topographical density variations, since lake beds represent large density contrasts between water and rock. However, our ability to investigate relations between lake bathymetry and gravity field effects is limited because of relatively sparse gravity data coverage over most lake surfaces, and the inability of satellite-derived gravity fields to capture the high frequency response to density effects. Satellite altimetry may provide an alternate method to determine a dense gravity field at surfaces of large inland lakes. We calculate satellite altimetry derived free-air (FA) anomalies over the Great Lakes in Canada, and evaluate them in terms of resolution and accuracy. The anomalies are compared to those calculated from the GRACE-derived EIGEN-GL04C gravity field to test for bias, and from shipborne or airborne results (where available) to test the accuracy of their determination of the high frequency component of the gravity field. Our results show that accuracies may reach less than 10 mgal, both in terms of bias and higher frequency effects, but are often tens of milligals. Furthermore, we find that altimetry does provide a higher resolution gravity field than satellite-derived gravity fields, and so if the accuracy can be improved in later efforts, it will be a useful tool for investigations of topographical density effects of lake waters.


IntroductIon
The study of Earth's gravity field, and the response of the gravity field to topographical density variations, has important implications in physical geodesy.One of these areas, and one not yet well explored, is the effect of topographical density variations on geodetic methods of height determination.Both geoid determination and orthometric heights depend directly on the behavior of Earth's gravity field, which in turn depends to a significant extent on the variations in crustal density near to the surface.While these effects have been approximated using a model of laterally-varying density (e.g., Huang et al. 2001;Hwang et al. 2003;Kingdon et al. 2005), the effects of density viewed as a three-dimensionally varying quantity are not well known.However, any study of the effect of these variations re-quires sufficient knowledge of the gravity field in a study area.One of the key study areas for investigations of threedimensional topographical density effects is the pronounced density contrasts existing at the bed of lakes.
Bathymetry can readily provide a three-dimensional model of these interfaces where available, but this is only useful if the gravity field around and over lakes is known to a sufficiently high resolution.While gravity anomalies can be computed from shipborne or airborne gravity data, in many areas this data is sparse.Also, while satellite gravity data, e.g., from GRACE, may furnish a long-wavelength field over lakes, high frequency density affects are indistinguishable in this signal, and so a higher resolution field is preferred if available.
Procedures for calculating gravity anomalies using satellite altimetry results and least squares collocation and Fourier transforms over oceans are well established (e.g., Sandwell and Smith 1997;Hwang et al. 2003;Sandwell and Smith 2005;Kim et al. 2008).These same procedures have been extended to inland water bodies (e.g., Hwang et al. 1998).Errors in either case result from: (1) corruption of waveforms by land interference and multi-returns from lake bottom for points in shallow waters near to shore; (2) errors in wet tropospheric corrections resulting from bad radiometer data; (3) error due to tidal effects; and (4) sea state bias.Over lakes, however, we expect (1) to have more influence than over oceans because lakewaters are overall more shallow.( 2) and (4) should be similar for both, while (3) should be slightly less in lakes since even in large lakes tidal influences are less than in oceans.Fortunately, the errors resulting from (1) are slightly diminished in the case of larger lakes, which are also the lakes most likely to have a significant effect on the gravity field.Furthermore, they can be mitigated in deep lakes by omitting satellite altimetry results near to shore (where waveforms are more erratic), and by including onshore gravity observations as part of the input to a least squares collocation procedure used to determine gravity anomalies (cf. Hwang et al. 2003).In some cases, such results have been extended to extrapolate information about density distributions over larger water bodies for some time now, in determinations of bathymetry from satellite altimetry results.Examples of such work are Smith and Sandwell (1994), Hwang (1999), and Vergos and Sideris (2002).Thus the relationship between density distributions and gravity, determined from satellite altimetry, is well established.What we hope to discover is whether such relationships can be investigated using satellite altimetry data for the topographical density distribution represented by lakewaters.
In this study, we attempt to calculate free-air gravity anomalies at a 1' 1' # resolution, better than that currently attainable from GRACE results alone, over the Great Lakes in Canada, with a view to finding sufficiently dense gravity anomalies for density investigations.The gravity anomalies are calculated using both Topex/Poiseidon and ERS-1/ Geosat derived height anomalies, in combination with freeair anomalies from terrestrial gravity data gathered from the United States Geological Survey's database in the United States, and the National Geological Survey's database in Canada.These are used as input to a least squares collocation, using the program colloc.fwritten by Prof. Chienway Hwang, and the output is examined by comparison with the GRACE-derived gravity field EIGEN-GL04C (Flechtner 2006), the EGM-96 field, and shipborne/airborne results to assess the accuracy of the gravity field determined from satellite altimetry.

Preparation of terrestrial Gravity data
Terrestrial gravity data for Canada and the United States is freely downloadable from the Canadian Geological Survey, and the National Geological Survey, respectively.Each agency provides files containing gravity point names, coordinates, dates and times of collection, raw gravity values, and standard deviations.Once complete data sets of raw gravity observations were downloaded, the data was filtered using MATLAB software to filter out data beyond 0.5 degrees of latitude and longitude of the Great Lakes, and those with standard deviations greater than 3 mgal.A visual inspection was used to verify that the data at least appeared to represent a continuous surface.
The remaining data points were converted into free air anomalies according to Hofmann- Wellenhof and Moritz [2005, Eqs.(3) -( 100)].Thereafter, the data set was exported from MATLAB and separated into free-air anomalies over lakes, to be used for later comparison with results; and terrestrial free-air anomalies, to be used as input for least squares collocation.
Finally, the terrestrial free-air anomalies were filtered using the box median filter included with the Generic Mapping Tools (GMT, Wessel and Smith 1999) to decimate it to a 1 arc-minute resolution.This filter divides input data into cells of user specified size, and removes from each cell all but the median data point.It was applied because our terrestrial gravity data was too dense, resulting in unnecessarily slow calculations and erratic behavior of our least squares collocation algorithm.Final data set coverage is shown in Fig. 1.

Preparation of ErS-1/GM and Geosat/GM data
ERS-1/GM and Geosat/GM data at a 5 Hz (1.4 km) sampling rate were retracked following the method of Sandwell and Smith (2005), to minimize the error in surface slope and to decouple it from wave height error.This is especially important because we use height gradients to derive gravity anomalies.Furthermore, Sandwell and Smith's method has proven useful in coastal areas, which are analogous to lake waters.
Because we expect altimetry results over lakes to contain some anomalous results not related to the mean lake surface, the retracked altimetry data was filtered for outliers by comparison of the observed heights with geoidal heights derived from the EIGEN-GL04C gravity field.Outlier testing was carried out using the sample standard deviation and sample mean for these differences, allowing for a bias between the altimetry results and the EIGEN-GL04C field.Pope's tau test was applied to create a confidence interval around the sample mean, and points where the difference fell outside of this interval were considered outliers.This filtering assumes that the surface of the lake follows a geopotential surface.Figure 2 shows data coverage before filtering, and Fig. 3 shows the outliers found.
It is evident from Fig. 2 that ERS-1/GM data alone does not offer better coverage than airborne or shipborne data.Geosat/GM data is necessary to improve upon the coverage from shipborne or airborne data.Also, a relatively small number of altimetry observations are removed as outliers, justifying our assumption that the lake surface coincides with a geopotential surface.
While many outliers are near to the shore, some occur near the middle of the lakes.Comparison of these with bathymetry reveals that they are often associated with formations on the lakebeds.These formations could influence the lake surface either by their interaction with water currents (which can be quite strong in the Great Lakes), or by the   way they effect the gravity field.In the latter case, they may represent interesting phenomena for topographic density effect investigations; however, without shipborne or airborne data-which is relatively sparse over the Great Lakes, as Fig. 1b shows-it is difficult to distinguish between these effects.
We can say that any local bathymetric formation would have to be quite large to have a noticeable effect on altimeter-derived gravity.For example, suppose a very optimistic sensitivity of altimeter results to gravity field variations of 3 mgal.This is suggested in Hwang et al. (2003) as the current best possible accuracy of altimetry-based gravity anomalies, based on results of e.g., Sandwell and Smith (1997), Andersen et al. (2001), andHwang et al. (2002).In order to cause an effect on the gravity field of this magnitude, a bathymetric formation would have to be very large.One possibility of many would be a 50-m high, 460-m wide, hemi-ellipsoidal formation of average crustal density (2670 kg m -3 ) at a depth of 100 m below the lake.Such a formation would be too small to have a distinguishable effect on a satellite-derived gravity field, such as the EIGEN-GL04C field used for outlier testing.However, it could affect the lake surface heights at up to two along track points, which might according to our testing be mistakenly considered outliers.
Finally, along track height gradients are calculated from the satellite altimetry-derived heights.Again assuming the lake surface to be an equipotential surface, these are shown to be the same as the geoid gradients, e.They are calculated according to: where ∆h is the height difference between two along-track points, and ∆d is the distance between the points.

calculation of Gravity Anomalies from terrestrial and Satellite Altimetry data
A remove-compute-restore technique is applied to calculate gravity anomalies.In this procedure, effects of a global geopotential model are removed from both the geoid gradients e and the terrestrial gravity anomalies ∆g.This is done according to: where C ge is the covariance matrix relating gravity to geoid gradients, C gg is the gravity-gravity covariance matrix, and is the geoid gradient-geoid gradient covariance matrix.The terms D g and D e represent noise in the gravity and geoid gradient observations.The output of the collocation could be at any suitable resolution, but in our case a resolution of 1' 1' # was chosen.The covariance functions needed in the least squares collocation computations are expanded into series of Legendre polynomials (Moritz 1980).The coefficients of the polynomials are derived from the error degree variances of the reference geopotential modes for degrees 2 to 360, and from Tscherning and Rapp (1974) Model 4 degree variances for degree 361 to infinity.The details are given in Tscherning and Rapp (1974).The actual computations of covariance functions used a modified version of a program in the computer package "GRAVSOFT".
After least squares collocation, the reference gravity anomaly values are restored to the final free-air anomaly values, to calculate final free-air anomalies: where ∆g ref is the EIGEN-GL04C-derived gravity anomaly, and the other quantities are as defined above.

Assessment of results
The accuracy of results was assessed by comparison with two global geopotential models, and with shipborne/ airborne gravity data.The shipborne and airborne data was converted into free-air anomalies using the same procedure as with the terrestrial data.Comparison with global models EGM96 and GGM02C was performed to search for biases in our results, and comparison with shipborne or airborne data was used to estimate errors in the high-frequency component of our results.Differences between the data sets were calculated, with the expectation that these differences would fall within the uncertainties of altimetry-derived gravity anomalies suggested by Hwang et al. (2003), in which case we could verify that the results are indistinguishable from each other.Also, the comparison with both high and low frequency data can give some idea where our data fits in the spectrum of resolutions of gravity data available over lakes.

comparison of Area results
The entire results of our gravity anomaly calculations are shown by the color plot in Fig. 4. For comparison, Great Lakes bathymetry retrieved from NOAA/GLERL (Schwab and Sellers 1996) is shown in Fig. 5.
At first glance, the results show a significant correlation with lake depths.The results over lakes range from -80 to 20 mgal, with the lowest values in shallow waters and the noisiest results in bumpy areas.This is an encouraging sign for the usefulness of these results for density investigations, because the effects of these density contrasts (between the lake bed and the water) are not otherwise available in these areas.Fields based on satellite observations, while they cover this area, are observed too far from the density contrasts to register their effects.The lack of response to highfrequency density contrasts is visible in the EIGEN-GL04C gravity field over the same area, expanded to degree 360, shown in Fig. 6.The EGM96 field is visually indistinguishable from the EIGEN-GL04C field at this resolution, and so is not displayed.
Although Figs. 4 and 6 show similar trends, the EI-GEN-GL04C field is much smoother than the satellite altimetry results.Also, while the EIGEN-GL04C model does show some response to bathymetry, high-frequency bathymetric phenomena are not distinguishable in the field.Thus satellite altimetry results provide additional high frequency information on topographic density formations neither available from low-resolution global geopotential fields, nor from low-frequency satellite gravimetry results.

comparisons Along test Profiles
Our results were next compared to free-air anomalies from shipborne or airborne gravity observations over the lakes.In this way, the high-frequency component of our results could be tested.Figure 7 shows the results for comparisons with all shipborne and airborne data over the lakes.The profiles are labeled P1 to P13, in the order in which they were tested.Profiles P2, P9, P10, and P5 are drawn from airborne gravity surveys, while the other profiles come from shipborne surveys.Although as can be seen above many of the airborne and shipborne results used for comparisons were within at most twenty mgal of the altimetry results, in some cases there were significant differences.These are usually in the case of shipborne results, resulting from a large bias in shipborne data as compared to both results from satellite altimetry and from geopotential models.Typical examples are profiles P1 for Lake Superior, and P4 for Lake Ontario, shown in Fig. 8. Gravity values given by EIGEN-GL04C are included in the plots (dotted line) alongside those from satellite altimetry (thick line) shipborne/airborne results (narrow line).This convention is continued in later profile plots.
These profiles show a distinct bias in the shipborne observations.Furthermore, P4 shows quite clearly the high frequency field visible from the satellite altimetry results.The wide dip near the right side of the graph is particularly telling, as it corresponds exactly with a dip in the bed of Lake Ontario.These results introduce an unexpected reason why satellite altimetry data is useful in studies of gravity over lakes: publicly available shipborne data is not all reliable.
In other areas, the shipborne data agreed well (within ~±10 mgal, without bias) with satellite altimetry results.This was usually the case with airborne data, such as along P2 and P9 in Lake Huron.In other cases, while airborne and altimetry results agreed very well with each other, they were different from EIGEN-GL04C results.This is the case for P10 over Lake Michigan.Both P9 and P10 are shown in Fig. 9 as illustrations of agreement between satellite altimetry and airborne results.
Note that the divergence between EIGEN-GL04C and airborne/altimetry results at the left side of Fig. 8b is likely the result of a long narrow dip in the bed of Lake Michigan which runs along the first half of the profile.Because of the narrowness of this formation, it does not appear in the EIGEN-GL04C field at only degree 360, although it has a noticeable effect on the other two sets of results.It is noteworthy that satellite altimetry and airborne results register this gravity disturbance with similar accuracy to each other.This implies a strong response to underwater density formations in satellite altimetry results.
Finally, there are some situations where the satellite altimetry results are simply inadequate, as revealed by biases from both the EIGEN-GL04C and shipborne gravity.Two such situations occur for profiles P6 and P13, shown in Fig. 10 below.
Sometimes there are clear physical explanations for these biases, as with P6 in Lake Erie, which shows a difference of about 30 mgal between satellite altimetry results and shipborne/EIGEN-GL04C results at its west side.In that case, the difference grows as the shallow lake (maximum depth of 64 m) becomes more shallow near to the shore.The difference is presumably a result of the effect of this shallowness on satellite altimetry results.Perhaps in the case of lakes-especially those with relatively shallow areas extending far from shore-it is better to decide which satellite altimetry data to use based on water depth, and not just proximity a shoreline.
In occasional situations, however, satellite altimetry data shows deviations from EIGEN-GL04C and shipborne results for no apparent reason.While this was rare among the profiles we tested, it is the case for P13, where the altimetry results differ by up to 22 mgal from the shipborne results and show a bias from EIGEN-GL04C results of about 15 mgal along most of the profile.For other profiles, all differences between sources of gravity data have been explicable for one of the reasons above: a bias in the airborne/shipborne data, the inability of the EIGEN-GL04C field expanded to degree 360 to register local density contrasts (those having high-frequency signals), and the errors in satellite altimetry results over shallow waters.None of these effects is clearly responsible for the differences between results visible along P13.The left half of the profile at least is over deep water, and so errors in altimetry results over shallow waters cannot have caused the bias.Comparison with bathymetry does not reveal any correlation with bathymetric formations, and even if it did we would expect to see these formations in shipborne results as well.And finally, the shipborne results agree well with EIGEN-GL04C, showing that it is the altimetry results which are biased.

Statistical Examination of results
Although our results were not rigorously tested, some comments can be made about their statistical behavior.In particular, the maximum, minimum, and average differences between our satellite-altimetry derived gravity field and both the EIGEN-GL04C and shipborne/airborne gravity results, and their standard deviations, can be examined.These statistics are provided in Table 1.
Notice first that the altimetry results are consistently lower than the EIGEN-GL04C gravity.This is anticipated, since the global geopotential model represents a smoothing of the gravity field, and thus is less affected by the density anomalies represented by the lakes.Because the water in a lake has a lower density than surrounding topography, it will reduce the value of gravity at the lake surface.However, while gravity at points on the lake's surface will be sensitive to these effects, a smoothed gravity field will be more strongly influenced by the more distant topographical masses surrounding lakes.Thus we expect the smooth EI-   GEN-GL04C field to be slightly stronger than one derived by satellite altimetry, seaborne, or airborne observations.If we follow the assumption that the EIGEN-GL04C field represents well the low-frequency gravity field over the lakes, we can consider the standard deviation of the differences between this and the altimetry results like an RMS error.If we do so, we can error in the differences of less than 13 mgal for most profiles, which is acceptable.However, the maximum and minimum values show that for most profiles there are some points which lie well outside of this range.
In constrast to differences with EIGEN-GL04C, differences with shipborne/airborne results do not show a variety of standard deviations.In some cases (profiles 4, 7, and 12), standard deviations of the differences are lower than for EIGEN-GL04C, suggesting that in these cases the shipborne/airborne data and the altimetry results both are modeling the high frequency component of the gravity field better than the EIGEN-GL04C model.In others (profiles 2, 3, 5, 8, 9, 10, and 13) the standard deviations are similar to those for the differences with the EIGEN-GL04C model.This is often the case for profiles based on airborne results (profiles 9, 10, and possibly 5), in which the airborne results tend to be smoother than the altimetry results.Finally, in some cases (profiles 1, 6, and 11) the standard deviations are much higher than those with the EIGEN-GL04C field.In these cases, it is often large peaks in the shipborne/airborne results which cause the difference.Note also that while some profiles show high magnitude minima for differences with EIGEN-GL04C, this is somewhat deceiving since most of these are the result of high frequency bathymetric variations rather than errors in the altimetry results.
Also, in general the altimetry results agree better with the EIGEN-GL04C model than with the shipborne and airborne results.This implies that there is little bias in the altimetry results (no more than 16 mgal in magnitude, for these profiles), while there is certainly a larger bias in some of the airborne/shipborne results (as much as 48 mgal).Profiles where bias is especially present in shipborne results are 1, 4, 5, 7, 8, and 11.However, for some profiles (2, 3, 9, and 10), there is a very good agreement with shipborne/ airborne results.For these profiles, the maximum and minimum differences are relatively small (they never exceed 20 mgal), while the standard deviation of the differences is even smaller (never exceeding 8 mgal) and the bias ranges from small (4 mgal) to nonexistent.It would be interesting to compare the altimetry results with the latest EGM model (EGM08, released 2008) over the Great Lakes, particularly the high frequency parts.

concLuSIonS
Satellite altimetry data provides a higher frequency gravity field over lake areas than satellite-derived geopotential models, and is sometimes better than publicly available shipborne gravity data in the same areas.Altimetry results are mostly unbiased with respect to satellite-derived global geopotential fields such as EIGEN-GL04C.Where a bias between the two is present, it is usual explicable by reference to bathymetry formations which are captured by the altimetry-derived gravity field, but not by the low frequency geopotential fields.Furthermore, the altimetry results agree often with airborne and sometimes with shipborne gravity results, showing that they are able to accurately model the high frequency component of the gravity field.In many cases where they do not agree with shipborne results, it appears that it is the shipborne data which is flawed, often having a constant bias from the EIGEN-GL04C field.However, the accuracy of altimetry results over lakes is on the order of tens of milligals, which is insufficient to distinguish many density formations.Noise of this order might easily be mistaken for lake density effects in the absence of data to validate the altimetry results.Furthermore, in some cases deviations of altimetry results from other data sources occur.Most of these are explicable due to shallow waters adversely affecting altimetry results, but in rare cases no explanation is apparent.Thus while altimetry promises to model well the high frequency gravity field in areas of sparse data over lakes, they still need to be refined before they can provide a reliable gravity field of high accuracy.

Fig. 4 .
Fig. 4. Results of free-air anomaly calculation for the Great Lakes.
In Eq. (2), for geoid gradients, e ref represents the gradients derived from a geopotential model-in our case, EIGEN-GL04C expanded to degree 360-and e res represents the residual geoid gradient.Likewise with Eq. (3), ∆g ref represents the EIGEN-GL04C gravity anomaly, and g res

Table 1 .
Statistics of differences in free-air anomalies between altimetry, EIGEN-GL04C, and shipborne/airborne results (all values in mgal).