WIP atmospheric and spherical cap corrections #168
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,13 +1,18 @@ | ||
| """ | ||
| Gravity corrections like Bouguer corrections, free air, spherical cap and | ||
| atmospheric. | ||
|
|
||
| Refer to Hinze(2005) for formula. | ||
| http://library.seg.org/doi/10.1190/1.1988183 | ||
|
|
||
| """ | ||
| import numpy as np | ||
|
|
||
| from .constants import GRAVITATIONAL_CONST | ||
|
|
||
|
|
||
| def bouguer_correction(topography, density_crust=2670, density_water=1040): | ||
| """ | ||
|
There was a problem hiding this comment. @craigmillernz any reason for this? The There was a problem hiding this comment. opps my mistake, I thought it a typo |
||
| Gravitational effect of topography using a Bouguer plate approximation | ||
|
|
||
| Used to remove the gravitational attraction of topography above the | ||
|
|
@@ -69,3 +74,72 @@ def bouguer_correction(topography, density_crust=2670, density_water=1040): | |
| density[oceans] = -1 * (density_water - density_crust) | ||
| bouguer = 1e5 * 2 * np.pi * GRAVITATIONAL_CONST * density * topography | ||
| return bouguer | ||
|
|
||
| def atmospheric_correction(topography): | ||
|
There was a problem hiding this comment. E302 expected 2 blank lines, found 1 |
||
| """ | ||
| Calculates atmospheric correction as per Hinze 2005 eqn 5. | ||
|
|
||
|
There was a problem hiding this comment. W293 blank line contains whitespace |
||
| .. math:: | ||
|
|
||
| g_{atm} = 0.874 - 9.9E-5 h + 3.56E-9 h^2 | ||
|
|
||
|
|
||
| in which :math:`h` is the elevation of the station above the ellipsoid | ||
| and :math:`g_{atm}` is the gravitational effect of the atmospheric the air | ||
| mass in mGal. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| topography : array or :class:`xarray.DataArray` | ||
| Topography height and bathymetry depth in meters. | ||
| Should be referenced to the ellipsoid (ie, geometric heights). | ||
|
|
||
|
There was a problem hiding this comment. W293 blank line contains whitespace |
||
|
|
||
| Returns | ||
| ------- | ||
| atmos_correction : array or :class:`xarray.DataArray` | ||
| The gravitational effect of the atmosphere in mGal | ||
|
|
||
| """ | ||
| atmos_correction = 0.874 - 9.9e-05 * topography + 3.56e-09 * topography**2 | ||
|
|
||
|
There was a problem hiding this comment. W293 blank line contains whitespace |
||
| return atmos_correction | ||
|
|
||
|
|
||
| def spherical_bouguer_cap_correction(topography): | ||
|
There was a problem hiding this comment. Does the density not enter into this equation? There was a problem hiding this comment. See comments above about using the full analytic expression. This approximation assumes 2670, so need to scale to other density. |
||
| """ | ||
| Calculates spherical cap correction for the Bouguer slab as per Hinze 2013. | ||
|
There was a problem hiding this comment. Is this a correction applied to the Bouguer correction or the effect of a spherical cap? In other words, does this require a Bouguer correction first or can be it directly applied to the gravity disturbance? There was a problem hiding this comment. It's my understanding that the spherical cap correction is applied to the Bouguer slab correction. The Bouguer correction has 3 parts, the infinite slab (2piGph, also known as Bullard A), the spherical cap correction (also known as Bullard B) and the terrain correction (Bullard C). However we can reformulate it to include the spherical cap correction in the Bouguer slab correction, so that there is only one "Bouguer correction" step to apply. |
||
|
|
||
|
There was a problem hiding this comment. W293 blank line contains whitespace |
||
| Gravity and Magnetic Exploration, Cambridge Press 2013 | ||
|
|
||
|
|
||
|
There was a problem hiding this comment. W293 blank line contains whitespace |
||
| .. math:: | ||
|
|
||
| g_{spher} = 0.001464139 h - 3.533047e-07 h^{2} + 1.002709e-13 h^{3} | ||
| + 3.002407E-18 h^{4} | ||
|
There was a problem hiding this comment. Could we use the analytical formula for the cap instead of this one? That way, the size of the cap can be a parameter instead of restricted to 167km (which might not apply for other planetary bodies or for satellite data). There was a problem hiding this comment. Yes that would be better, you mean as defined by LaFehr 1991 (Geophysics 56)? That allows for density, cap radius and earth radius to be specified.
The approximation assumes a constant density of 2670 kg/m3, so when a different reduction density is used the affect is more apparent. shown here for 2400 kg/m3 density when allowing for different ellipsoids and cap extents. There was a problem hiding this comment. reformulated to include the spherical cap correction in the "Bouguer slab formulation" to make it one step The spherical_bouguer_cap_correction now returns both the correction and the spherical corrected infinite slab and the logic check between the infinite slab value (55.971) plus the spherical cap correction (0.644) equals the spherical bouguer correction (56.615) checks out. There was a problem hiding this comment. Let me know your preference for a 1 step or 2 step correction. There was a problem hiding this comment. Right, I makes more sense now. I'm leaning towards a 1 step process since this is only ever used with the Bouguer correction anyway. We could include an argument like There was a problem hiding this comment. I agree. What I proposed is a single function called |
||
|
|
||
| in which :math:`h` is the elevation of the station above the ellipsoid | ||
|
There was a problem hiding this comment. This strikes me as a bit odd since the effect of a cap depends on 2 things: the height of the observation point and the thickness of the cap (i.e., height of topography). Is this equation assuming that they are the same? In fact, it should actually depend on the radial coordinate for both of these, not really the heights (which the authors probably assumed as some mean Earth radius). This will cause the equation to fail for lunar or martian applications, where the reference radius is much smaller. This can quickly get very complicate if we assume that the heights are relative to the ellipsoid not a sphere. One compromise could be:
There was a problem hiding this comment. See comments above about using the full analytic expression. I think if we use that then we can include both the planet radius and cap extent (though there is pretty good evidence that beyond 167 km it gets more imprecise as average crustal densities are likely to vary over large scales.) See Figure 4 in LaFehr 1991. The full expression does assume a constant radius and is most valid at mid latitudes where radius changes are small. If we have a method of calculating the specific raidus at the station location, then that could be used instead, but Figure 6 from LaFehr shows that it's only at high altitudes (4000m +) that the effect becomes more than instrument error (i.e max of ~20 uGal difference between pole and equator at 4000m), so unless you have a survey that covers the hemisphere I don't think it matter. There was a problem hiding this comment. I'm still a bit confused here about the radius referred in the papers. Is that assuming that observations are located at the surface of topography? In that case, the radius of the observation point and the radius for the spherical topography would be the same. But if we have observations from airplanes or satellites, this wouldn't work. Right? But then again, the Bouguer correction itself doesn't work in these cases anyway. Either way, I like the full analytic expression! It's good to have these things be explicit. I just had a bug in some code doing a Bouguer correction for Mars and forgetting to set There was a problem hiding this comment. I'm not too familiar with airborne data processing (and even less satellite) but i did some googling and it seems airborne gravity is processed similar to ground gravity (except it requires the EotVos correction of course, which we could include as a separate function). See report below for processing steps: http://www.geosciencebc.com/i/project_data/QUESTdata/GBCReport2008-8/Gravity_Technical_Report.pdf There is a "curvature" correction that requires the input of the ground surface (not the aircraft altitude). When I computed this it is approximately the same as the spherical_bouguer_cap_correction, (0.64368 mGal). So i think that for airborne data you would calculate the "normal gravity" at the aircraft elevation (which is where the sensor is) and then calculate the spherical bouguer slab at the elevation of the ground below the aircraft (which can be dervied from the laser altimeter data etc (flight height - laser altimeter = ground height). In summary, it seems we can use the full spherical bouguer for both ground and airborne gravity. The user needs to ensure they use the correct elevations (aircraft vs ground) for the various parts of the complete Bouguer Anomaly formulation. There was a problem hiding this comment.
Yes. The problem I mentioned is that gravity should decay with observation height but in all of these different Bouguer corrections the observation height isn't a parameter. That is because the Bouguer plate has constant gravity with height due to being infinite. This is an OK approximation if the observation is on top of the topography, but it's an over-correction when the observations are at greater height (airborne). I was wondering if the spherical version took this into account. But apparently it doesn't.
Yes! This is why we're avoiding the word "height" when talking about Bouguer and using "topography" instead to make it clear that this is not the observation height but the topographic height. |
||
| and :math:`g_{spher}` is the gravitational effect of the spherical cap to | ||
|
|
||
| 166.7 km in mGal. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| topography : array or :class:`xarray.DataArray` | ||
| Topography height and bathymetry depth in meters. | ||
| Should be referenced to the ellipsoid (ie, geometric heights). | ||
|
|
||
|
There was a problem hiding this comment. W293 blank line contains whitespace |
||
|
|
||
| Returns | ||
| ------- | ||
| spher_cap_corr : array or :class:`xarray.DataArray` | ||
| The gravitational effect of the spherical cap in mGal | ||
|
|
||
| """ | ||
| spher_cap_corr = 0.001464139 * topography - 3.533047e-07 * topography**2 | ||
|
There was a problem hiding this comment. E222 multiple spaces after operator |
||
| + 1.002709e-13 * topography**3 + 3.002407E-18 * topography**4 | ||
|
|
||
|
There was a problem hiding this comment. W293 blank line contains whitespace |
||
| return spher_cap_corr | ||
|
|
||
|
|
||
|
|
||
|
|
||
|
There was a problem hiding this comment. W293 blank line contains whitespace |
||
|
|
||
|
There was a problem hiding this comment. W293 blank line contains whitespace |
||
|
|
||
|
There was a problem hiding this comment. W391 blank line at end of file |
||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
W291 trailing whitespace