Determination of the Topography-Bounded Atmospheric Gravity Correction for the Area of Poland

The atmospheric correction (Eq. 1) is the difference of two components. The first is the normal gravity component defined by the atmospheric masses contained within the reference ellipsoid (the normal gravity component g^NA), which can be defined as (e.g. Sjöberg, 1993): (2) gNA=GMar2 {g^{NA}} = {{G{M_a}} \over {{r^2}}} where G is Newton’s gravitational constant, M_a is the mean mass of the Earth’s atmosphere and r is the geocentric radius of the computation point.

The second component (spherical atmospheric component g^SA) is determined by the atmospheric masses constituting spherical, homogeneous layers extending from the surface of the reference sphere as the lower boundary, and it can be written as (3) gSA=GMr−Rr2 {g^{SA}} = {{G{M_{\left( {r - R} \right)}}} \over {{r^2}}} where M_(r−R) denotes the sum of the masses of homogeneous, spherical layers of the atmosphere located between the reference sphere of the radius R and the computation point.

Hence, in general, we will define Eq. (1) as (4) δgatm=gNA−gSA \delta {g_{{\rm{atm}}}} = {g^{NA}} - {g^{SA}}

If topography is adopted as the lower boundary of the atmosphere, the gravity effect of atmospheric masses (g^ETA) will be defined as (Figure 1) (5) gETA=−∂VETA∂r {g^{{\rm{ETA}}}} = - {{\partial {V^{{\rm{ETA}}}}} \over {\partial r}} where (6) VETA=G∭ΩρAldVΩ {V^{{\rm{ETA}}}} = G\mathop{\int\!\!\!\int\!\!\!\int}\limits_{\kern-5.5pt \Omega } {{{{\rho _A}} \over l}d{V_\Omega }}

In Eq. (6), Ω is the volume of integration, ρ_A is the atmosphere density distribution function, l is the distance between the attracting masses and the attracted point and dV_Ω is the element of volume. The Ω volume covers the masses of the atmosphere from the surface of the topography to a height at least equal to the highest mountain peaks.

Alternatively, the g^ETA component can be determined from the difference (Mikuška et al., 2008) (7) gETA=gSA−gTA {g^{{\rm{ETA}}}} = {g^{SA}} - {g^{TA}} where g^TA is the topographic–atmospheric correction equal to the gravity of all Earth topography, assuming its density is equal to the density of the atmospheric masses (Figure 1): (8) gTA=−∂VTA∂r {g^{TA}} = - {{\partial {V^{TA}}} \over {\partial r}} where (9) VTA=G∭ΘρAldVΘ {V^{TA}} = G\mathop{\int\!\!\!\int\!\!\!\int}\limits_{\kern-5.5pt \Theta } {{{{\rho _A}} \over l}d{V_\Theta }}

In Eq. (9), Θ represents the volume of all Earth topography.

Figure 1.

The Ω and Θ volumes adopted to determine the g^ETA component

Now, the topography-bounded atmospheric correction will be defined as (10) δgatmETA=gNA−gETA \delta g_{{\rm{atm}}}^{{\rm{ETA}}} = {g^{NA}} - {g^{{\rm{ETA}}}}

Let us note that taking into account the computational costs, determining the g^ETA component according to Eq. (7) is more advantageous than according to Eq. (5) due to the much smaller integration domain (in areas of seas and oceans, there is no topography, so they can be omitted). For this reason, we used it in further investigations.

The method of determination of Eq. (8) depends on the assumptions made regarding the coordinate system, available digital elevation model (DEM) and the atmospheric density distribution model. Assuming a spherical, geocentric coordinate system, DEM will be defined as a regular grid of spherical tesseroids.

By assigning each tesseroid a constant density ρAi \rho _A^i , the g^TA component can be determined based on (11) gTA=∑i=1nδgi {g^{TA}} = \sum\limits_{i = 1}^n {\delta {g_i}}

In Eq. (11), n is the number of tesseroids used in the calculations and δg_i is the component of gravity determined based on a single tesseroid defined as (Heck and Seitz 2007) (12) δgi=GρAi∫λ1iλ2i∫φ1iφ2i∫r1ir2ir−r′cosψr′2cosφ′l3dr′dφ′dλ′ \delta {g_i} = G\rho _A^i\int\limits_{\lambda _1^i}^{\lambda _2^i} {\int\limits_{\varphi _1^i}^{\varphi _2^i} {\int\limits_{r_1^i}^{r_2^i} {{{\left( {r - r'\cos \psi } \right)r{'^2}\cos \varphi '} \over {{l^3}}}dr'd\varphi 'd\lambda '} } } where φ′, λ′ and r′ are coordinates of the running integration point; λ1i \lambda _1^i , λ2i \lambda _2^i , φ1i \varphi _1^i , φ2i \varphi _2^i , r1i r_1^i and r2i r_2^i are the coordinates defining thesseroid i of DEM and ψ is the angle between the computation point (with coordinates φ, λ, r) and the running point: (13) cosψ=sinφsinφ′+cosφcosφ′cosλ′−λ \cos \psi = \sin \varphi \sin \varphi ' + \cos \varphi \cos \varphi '\cos \left( {\lambda ' - \lambda } \right)

Because there is no solution of the integral in Eq. (12) and the area of integration is very large and includes very distant masses, the δg_i components were determined using three approximations. In the immediate vicinity of the calculation point (up to 5 km), tesseroids were replaced by rectangular prisms of the same mass, for which a closed solution of the appropriate integral is known (Nagy et al., 2000). For tesseroids located further away (5 − 500 km), we used the approximate solution of the integral in Eq. (12) in the form proposed by Heck and Seitz (2007). For further masses, tesseroids were replaced by point masses located in their centres.

Due to changes of atmospheric density in the vertical direction, determination of the g^TA component will also require vertically dividing individual DEM blocks into segments of appropriate height.

The United States Standard Atmosphere 1976 model (USSA76) was used to estimate the density of the atmosphere, ρAi \rho _A^i , needed to implement Eq. (12). Based on the tabular values included in model USSA76, an analytical model was determined using the least squares method in the form: (14) ρAH=∑i=04aiHi {\rho _A}\left( H \right) = \sum\limits_{i = 0}^4 {{a_i}{H^i}} where ρ_A is determined in kg/m³, H is the height in metres and a_i is a coefficient given in Table 1.

Based on the model (Eq. 14) and assuming H = r − R, the M_(r−R) value can be written as the integral: (15) Mr−RH=4π∫0HR+H2ρAHdH=4π∫0HR+H2∑i=04aiHidH {M_{\left( {r - R} \right)}}\left( H \right) = 4\pi \int_0^H {{{\left( {R + H} \right)}^2}{\rho _A}\left( H \right)dH} = 4\pi \int_0^H {{{\left( {R + H} \right)}^2}\left( {\sum\limits_{i = 0}^4 {{a_i}{H^i}} } \right)dH}

Its solution is given in the form: (16) Mr−RH=4π∑i=04si {M_{\left( {r - R} \right)}}\left( H \right) = 4\pi \left( {\sum\limits_{i = 0}^4 {{s_i}} } \right) where:

• s0=a0HR2+RH+13H2 {s_0} = {a_0}H\left( {{R^2} + RH + {1 \over 3}{H^2}} \right) ;

• s1=a1H212R2+23RH+14H2 {s_1} = {a_1}{H^2}\left( {{1 \over 2}{R^2} + {2 \over 3}RH + {1 \over 4}{H^2}} \right) ;

• s2=a2H313R2+12RH+15H2 {s_2} = {a_2}{H^3}\left( {{1 \over 3}{R^2} + {1 \over 2}RH + {1 \over 5}{H^2}} \right) ;

• s3=a3H414R2+25RH+16H2 {s_3} = {a_3}{H^4}\left( {{1 \over 4}{R^2} + {2 \over 5}RH + {1 \over 6}{H^2}} \right) ;

• s4=a4H515R2+13RH+17H2 {s_4} = {a_4}{H^5}\left( {{1 \over 5}{R^2} + {1 \over 3}RH + {1 \over 7}{H^2}} \right) .

Eq. (16) allows determining the exact values of the g^SA component (Eq. 3), taking into account vertical atmospheric density changes in form of the model given by Eq. (14).

As delineated in the preceding section, owing to variations in atmospheric density along the vertical, the computation of the g^TA component according to Eqs (11) and (12), also requires the vertical partitioning of each individual DEM block into segments of suitable heights and constant density. To assess the vertical resolution of such segmentation, we will use the g^SA component, whose exact values can be determined. Let us note that the g^SA component (Eq. 3) can also be represented in an approximate form through the division of the atmosphere into m homogenous layers characterised by constant density and altitude. The M_(r−R) value can, therefore, be expressed in the following form: (17) Mr−Rr≅43π∑j=1m−1ρAjrj+13−rj3 {M_{\left( {r - R} \right)}}\left( r \right) \cong {4 \over 3}\pi \sum\limits_{j = 1}^{m - 1} {\rho _A^j\left( {r_{j + 1}^3 - r_j^3} \right)} where r_j = R + dh × (j − 1), dh=r−Rm dh = {{r - R} \over m} is the height of a single layer and ρAj \rho _A^j is the density of the atmosphere in the middle of the height of the layer j determined on the basis of Eq. (14).

The error in determining the g^SA component caused by using the approximate M_(r−R) value defined by Eq. (17) instead of the exact value provided by Eq. (16) can be estimated by considering the difference: (18) δgSA=gapproxSA−gexactSA \delta {g^{SA}} = g_{{\rm{approx}}}^{SA} - g_{{\rm{exact}}}^{SA} where gapproxSA g_{{\rm{approx}}}^{SA} and gexactSA g_{{\rm{exact}}}^{SA} represent the g^SA components determined using the M_(r−R) values in their approximate form (Eq. 17) and exact form (Eq. 16), respectively.

To estimate the δg^SA errors for different dh, the gapproxSA g_{{\rm{approx}}}^{SA} values were determined using Eqs (3) and (17) for one point situated at a height of 2,500 m (the highest point in the area of Poland) and for another point situated at a height of 8,500 m (the highest point on the Earth). Similarly, using Eqs. (3) and (16), the gexactSA g_{exact}^{SA} values were determined for the same points. The determined δg^SA values are presented in Figure 2.

Figure 2.

The δg^SA values for points with H_p = 2,500 m (red line) and H_p = 8,500 m (black line). The symbols in both charts indicate the dh values for which the calculations were performed.

Initially, it is noteworthy that the δg^SA error values are almost equal to 0 for dh ≤ 100 m (in both analysed cases), whereas for dh > 100 m, they are noticeably larger for the point located at a higher elevation. This indicates that for points with heights lower than those examined, the δg^SA errors will not exceed the levels depicted in Figure 2. For instance, for points situated within the geographic boundaries of Poland (not exceeding 2,500 m in altitude), assuming dh = 500 m, the gapproxSA g_{{\rm{approx}}}^{SA} component can be computed with the δg^SA errors not exceeding 0.02 μGal (as depicted by the red graph). Whereas, for locations not exceeding 8,500 m in elevation (i.e. points on the Earth's surface), the δg^SA errors may ascend to 0.05 μGal (as indicated by the black graph). Given that our calculations maintain an accuracy of approximately 1 μGal, such discrepancies remain within acceptable bounds.

The above discussed g^SA component errors also mirror those of the g^TA component, determined with the same vertical division of tesseroids constituting DEM. These errors will be equal if the topography consists of layers with constant heights (either 2,500 or 8,500 m, respectively). Consequently, the ascertained δg^SA errors for a height of 8,500 m can be considered the maximum errors in determining the g^TA component (for any point on the Earth's surface), resulting from the need to take into account changes in atmospheric density in the vertical direction by dividing tesseroids that constitute DEM into segments with constant of both atmospheric density and dh values. Considering the topography in the study area is generally lower overall, with the highest mountains situated at considerable distances, covering only a fraction of the Earth's surface, the error resulting from adopting dh = 500 m for the vertical tesseroids division will be substantially lower than the 0.05 μGal pointed above. Hence, dh = 500 m was selected to compute the g^TA component. Naturally, the g^SA component was determined in its exact version, based on Eqs (3) and (16).

The calculations were based on two DEMs: the SRTM v4.1 (3″ × 3″ model) (Jarvis et al. 2008) and the ETOPO1 (1′ × 1′ model) (NOAA National Geophysical Data Center 2009). Based on these models, five sets of DEM were prepared and directly used in the calculations for various distances from the calculation point:

• 3″ × 3″ - model for calculations at a distance of up to approximately 10 km;

• 15″ × 15″- model for calculations at a distance of approximately 10 − 50 km;

• 30″ × 30″- model for calculations at a distance of approximately 50 − 167 km;

• 1′ × 1′- model for calculations at a distance of approximately 167 − 500 km;

• 6′ × 6′- model for calculations carried out for topographic masses located further than 500 km.

The primary study area was located between parallels 48.4° − 55.1° and meridians 13.7° − 24.5° (Figure 3a). For this area, the calculations were performed for a regular grid of points with a resolution of 93″ × 90″. The area is mostly lowland. In the north, it covers a part of the Baltic Sea and in the south, it covers the Carpathian Mountains and the Sudetes (in the southwest). The highest region in the analysed area is the Tatra Mountains, with peaks reaching 2,499 m (Rysy) for the Polish part and 2,655 m (Gerlach) for the Slovak part of the Tatra Mountains. For part of the Tatra Mountains, marked in Figure 3a with a red rectangle, calculations were performed in a grid of points with a resolution of 12″ × 9″. This area is shown in Figure 3b.

Figure 3.

Relief maps of the study areas for a) the whole areas of elaboration and b) the Tatra mountains