Geology

The generalized geology, modified from Wager and Deer (1947), shown in figure 1 consists of a high-grade Archean gneissic metamorphic terrane (unit G n) partially overlain by Cretaceous sediments (unit S) and Tertiary basaltic flows, tuffs, sills, and dikes (unit B a s). The Skaergaard intruded this terrane during the Eocene together with some later sills and dikes. The entire section has been tilted in a broad flexure to the southeast from about 10° in the northern part of the intrusion to 30° in the Tertiary volcanic rocks south of the intrusion (Wager and Brown, 1967). Restoration of tilting indicates that the original zonal layering of the intrusion is essentially bowl shaped (Gettings, 1976). The intrusion is mapped into upper, middle, and lower zones identified by crystal composition (fig. 2; Gettings, 1976).

In addition to the layered rocks grouped into the lower, middle, and upper zones that make up the bulk of the intrusion, a marginal border group including the chilled margin that is up to several hundred meters thick and an upper border group are mapped (fig. 2). On the east side a younger, thick gabbro dike, named the Macrodyke (M, fig. 2), intrudes a short distance into the east side of the Skaergaard intrusion but is believed to continue beneath part or all of the Skaergaard (McBirney, 1996). A picrite body older than the Skaergaard intrusion is mapped to the north of the intrusion (P on fig. 2), as well as several older and younger gabbro dikes and sills in the surrounding terrane are not a part of the Skaergaard intrusion (McBirney, 1996). Subsequent to its solidification, a number of thin basaltic dikes intruded the Skaergaard gabbro and the surrounding host rocks. These dikes are related to the flexuring of the area (Wager and Brown, 1967). Figure 3 shows a cross section (profile A–A′) from north to south across the intrusion. The dolerite Basistoppen sheet intrudes the southern part of the Skaergaard intrusion beneath the upper border group (fig. 3) (Wager and Brown, 1967).

The data collection and processing followed by the results of the modeling of the gravity and aeromagnetic data are described below. Both two-dimensional (2D) and three-dimensional (3D) modeling of the gravity anomaly map and 3D modeling of the aeromagnetic data were carried out using models that incorporated the dipping layering of the intrusion. These models required derivation and coding of computer programs for cylinders and polygonal layers that are dipping relative to horizontal. All computations except for the presentation of the gravity and aeromagnetic anomaly maps as shaded relief images were carried out with an IBM 360/50 at the University of Oregon Computing Center during the period 1971–1973. Computer input at that time was by decks of punched cards, and output was punched cards and printed paper. Only some of the computer printouts were extant at the time of the writing of this report. These datasets have not been superseded by publicly available modern data and are available in the accompanying data release by Gettings and Parks (2025).

Gravity Survey

The Skaergaard gravity survey field work was completed in the period from July 24 to August 16, 1971, by H.R. Blank and M.E. Gettings using LaCoste-Romberg G-series gravimeter G-141 and Worden Master gravimeter W533. A total of 168 stations were established, of which 82 had altitudes determined by altimetry and 86 had known altitudes, mainly at approximate sea level determined by visual observation on the shoreline. Horizontal locations were recorded by labeled pinpricks on the 9″×9″ stereo aerial photographs. One day of helicopter use enabled us to obtain stations on several of the high, rugged basalt peaks east and south of the intrusion and Gabbrofjaeld (peak 1360, fig. 1). Approximately 14 traverses by foot, several including fly camps for two or three days were carried out, and at least 9 boat traverses were done. Several boat traverses were carried out with a motorized lifeboat and a rowboat. The expedition ship, the Signalhorne, was used for transport to and from various foot traverses and across Kangerdlugssuaq Fjord to Amdrup Point and Kap Deichmann. Watkins Fjord was obstructed by ice during the expedition and no stations could be established on the north shore. The captain of the Signalhorne plotted depth soundings on a map throughout the expedition to provide data on water depths along the coastlines, Kangerdlugssuaq and Mikis Fjords, and in Uttental Sund. Gravity station locations were compiled on a 1:50,000-scale base map (described herein) by scaling from features that could be located on both the aerial photographs and the base map. The only topographic map available at the time was the 1:250,000-scale map sheet Gronland 68.0.3 Kangerdlugssuaq of the Danmark Geodaetisk Institut (1969). This map was made from trimetrogon aerial photography taken by the United States Army Map Service (AMS) in the period immediately after World War II, coupled with ground control from surveys in 1933–34 and place names from 1955 expeditions. The base map was made from a photographic enlargement of the map sheet to 1:50,000-scale for an area 43 km east-west by 38 km north-south centered on the Skaergaard intrusion. The AMS trimetrogon aerial photographs, consisting of a right and left side-looking oblique photograph and a vertical photograph with stereo overlap, were obtained and used in the field for station locations. The resulting station distribution is 1 gravity station per 4 square kilometers (km²).

The gravity survey data were processed during 1971–1972 by M.E. Gettings and H.R. Blank, with some help from students in the digitization of 1:250,000-scale topographic maps for terrain corrections. It should be noted that the observed gravity values are on an arbitrary datum based on the reading of LaCoste-Romberg gravimeter G-141 on a base station upon return to the United States with no drift control. Drift of both gravimeters was computed using assumed linear drift between closures on the Skaergaard base station or temporary bases such as when working in the Mikis Fjord–Sodalen Valley area. Altitudes were determined by altimetry except those stations at or near sea level that were determined by Brunton compass leveling. Considerable effort was put into processing the altimeter data to obtain the best possible estimates of altitude. Given large barometric pressure variations due to changing weather conditions, most altimeter stations employed two or three altimeters so that a better estimate could be obtained. The correction scheme outlined in appendix 1 was used. Additionally, the barometric pressure records at three-hour intervals from the Danish weather station on the island of Nordre Aputateq, about 45 km south-southwest of the Skaergaard area, were obtained, which enabled reasonably good barometric corrections. After all corrections, the distribution of maximum difference between altimeters at the stations is shown in figure 4. The quantiles for these differences were (in feet [ft]) 0 percent, 0; 25 percent, 9; 50 percent, 18; 75 percent, 27; and 100 percent, 79, suggesting that about half of the altimetric altitudes are accurate to about ±10 ft.

Terrain corrections were carried out in four parts. The available 1:250,000-scale topographic maps were digitized using a grid of 1′ latitude by 2′ longitude for the survey area and at 3′ by 6′ intervals out to 167 km beyond the survey area with the estimated centroid altitude. Because bathymetric data were not available, the digitization was done at the ice altitude and sea surface. The Hayford-Bowie scheme of zone radii was used (Swick, 1942; Hayford and Bowie, 1912). The outer zone corrections from 2.29 km to 21.0 km were done using the 1′×2′ digitization, and corrections from 21.0 km to 166.7 km using the 3′×6′ digitization were computed using the program by Plouff (1977). The corrections for all stations from 68 m to 2.29 km were computed by transparent template over the 50,000-scale topographic map by H.R. Blank. The innermost corrections from the station to a radius of 68 meters (m) were computed by M.E. Gettings from the field notes using a combination of Hayford-Bowie quadrants and displaced half-slopes. Finally, the 1′×2′ prisms were given average density values based on the geologic map (Wager and Brown, 1967) and the measured densities of hand specimens described below, and corrections to the usual constant density terrain corrections were computed. Figure 5 shows the frequency distribution for the final terrain corrections. The quantiles for the terrain corrections were (in milligals [mGal]) 0 percent, -1.3; 25 percent, 3.2; 50 percent, 5.7; 75 percent, 8.9; and 100 percent, 37.5. The lack of topographic maps at scales larger than 1:250,000 meant that the inner zone terrain corrections, especially the Hayford-Bowie A–F zone (0–2.29 km about the station), have a large possible error. Given this, the error envelope for the terrain corrections is probably about ±20 percent. Thus, the error for half of the stations is about ±1.2 mGal or less, with larger uncertainties for stations with larger corrections.

As previously noted, the terrain corrections were carried out with data digitized to the ice altitude and sea surface because there was no bathymetric or ice thickness data. To evaluate the uncertainty introduced by this, several computational experiments were carried out. A cross section consisting of a portion of an ellipse terminated by a horizontal plane (see appendix 2) was chosen as a model of glacial valleys (fig. 6) and was used to estimate glacier thicknesses and fjord depths. Figure 7 compares the elliptic fit with soundings taken by the expedition ship on a traverse across Kangerdlugssuaq Fjord. For glacier profile estimates, the seaward part of young glacial Sodalen Valley (location in fig. 1) is bedrock, and a topographic profile across it from the base map was used to define the parameters.

Corbató (1964) derived the expression for the gravity effect of a 2D model using a cross section with a vertical axis parabola for the base truncated by a horizontal line for the top, and Corbató (1965) applied it to a glacier in Washington State where drill hole depths were available. M.E. Gettings compared this model to the similar model of a cross section with a portion of an ellipse and concluded the ellipse was a better representation as it gives a better fit to the increasing curvature at the model edges, even though the ellipse model case requires numerical integration. Using the width-to-depth ratios and slope at the edge from the Sodalen Valley profile (0.125 and 0.465), the ellipse semi-major and semi-minor axes were computed from the formulae in appendix 2, and then the gravity effect was computed from the formulae in appendix 3 for a range of glacier widths from 0.25 to 5.00 km. The results are summarized in figure 8 for a gravity station on bedrock at the edge of the glacier and for a station 100 m away from the edge on bedrock. Almost all the glaciers with stations near their edge in the survey were less than a kilometer in width, so the error in not correcting for the ice is between −0.2 and −2 mGal. The fjords have a greater effect as they are generally much wider; however, the big fjords are 5 or more km distant from the intrusion except on the southwest and had an approximately constant effect that did not greatly affect the residual anomaly due to the intrusion. On the southwest, the water is fairly shallow (reaching 200 m depth 1–2 km offshore), and in Uttental Sund along the west margin of the intrusion, maximum depths are 50 m or less except for a narrow 100 m depth in the southern part opposite Forbindelses Glacier (Forbindelsesgletscher), so the error in the gravity anomaly is not large and is almost systematic because of its location along the west side of the intrusion. The effect of correcting for the water areas would be to increase the residual anomaly somewhat because the maximum residual anomaly (discussed in the “Gravity Survey Analysis” section) falls in part along southern Uttental Sund and Skaergaardsbugt (the peninsula south of Uttental Sund, fig. 2). The maximum effect from the water is estimated to be 1.5 mGal.

Standard Bouguer gravity anomaly reduction formulae were used to compute the Bouguer gravity anomaly at reduction densities of 2.67 and 2.77 grams per cubic centimeter (g/cm³). Although the median density for the measured hand specimens of gneiss is 2.675 g/cm³ (refer to table 1) a higher density of 2.77 g/cm³ was used because the gneiss contains numerous mafic dikes and gabbro bodies to the north and west of the Skaergaard intrusion. East and south of the intrusion it is bordered by gneiss overlain by basalt that has a median density of 2.95 g/cm³. The resulting anomaly map was still strongly influenced by the dense basalts and gabbros to the east and south, and because an interpretation objective was to use the software of Cordell and Henderson (1968) from sea level downward, it was decided to model all terrain above sea level in the gravity reduction. This is known as a variable density Bouguer gravity reduction (Vajk, 1956; Grant and Elsaharty, 1962; Gimlett, 1964).

The variable density reduction was accomplished in several steps. A program allowing different densities for each prism was coded similar to that of Plouff (1977) for the inner zone corrections (out to 15 km about the station) using the vertical line element prism formula outside of 2.29 km from the station but with density specified for each prism and a circular cylinder with specified density inside the 2.29 km radius. Using the geologic map (Wager and Brown, 1967) to determine rock unit, density data were measured from the 238 hand specimens and the 38 published in Wager and Deer (1939), and the 200-m topographic contours, the equal mass elevation and density contrast from 2.77 g/cm³ was estimated for each of the 1′ latitude by 2′ longitude prisms (fig. 9) within 15 km of stations. The hand terrain corrections (0–2.29 km about the station) were multiplied by the ratio of the density estimated at the station divided by 2.77 g/cm³. Repeating the hand terrain corrections on a variable density map was investigated for two stations in rugged basalt topography, and the difference from the ratio correction was in the hundredths mGal, so it was decided not to recalculate the hand corrections with the variable density template. The resulting gravity anomaly values thus have been corrected for the geology above sea level and reflect density variation from sources below sea level. The final corrected station gravity anomaly values were used to generate a grid using 1/r² interpolation factor for surrounding stations where r is the distance from the grid intersection to the gravity station. The grid had a 1 km × 1 km cell size with an arbitrarily chosen origin 12.20 km south and 7.05 km west of peak 1277 (fig. 1). Figure 10 shows the final variable density Bouguer gravity anomaly map as a shaded relief image. This image used a new grid made using minimum curvature gridding techniques (Briggs, 1974) and was made using Oasis Montaj software (Geosoft, 2015). The variable density gravity anomaly map of the Skaergaard and its local surrounding area was originally published in Gettings (1976).

Aeromagnetic Survey

The aeromagnetic survey was conducted on July 5–7, 1971, using a Varian proton precession magnetometer with sensor towed by a Beechcraft Twin Bonanza aircraft when the weather cleared for a 2.5-day period. Flights were conducted from Isafjordur, Iceland, the northwesternmost airstrip available with logistic support that was 450 km from the study area. The flight on July 5, 1971, was flown at 15,000 ft altitude to take vertical photographs using a Hasselblad camera that was used to plot and locate the tracking photographs from the magnetometer flights the following 2 days (changes in the glacier shapes and snow cover from the stereo photography meant that the low level photographs from the magnetometer flight could not be reliably located on the stereo photographs without first locating the higher altitude photographs on the stereo photographs). The magnetometer flights were flown at a barometric altitude of 5,000 ft on July 6 and 7, 1971. Nominal flight line spacing was 1 km and navigation was visual. Four north-south tie lines were flown, one near the east and one near the west edge of the survey area, and two were flown in the center about 2 km apart, both crossing the intrusion. The main survey was flown in twenty east-west flight lines, with two approximately 10-km north-south lines over the northern two thirds of the intrusion north to the Watkins Fjord shoreline. Tracking photographs were taken with 35mm cameras vertically through a port in the aircraft belly and magnetometer data were recorded on a strip chart recorder with time and navigation fiducial marks.

Diurnal variation was estimated using the records of the Leirvogur Magnetic Observatory in Iceland (lat 64.183° N., long 21.70° W.), the nearest station, 620 km to the southeast. Vertical tracking photographs from the aircraft were plotted on the Hasselblad photographs and then using the Hasselblad photographs the points were located and plotted on the 9″×9″ stereo photographs. Using the stereo photographs, the points were then compiled on a base map at 1:50,000 scale. Magnetometer readings from the magnetometer strip chart were measured and plotted at approximately 300- to -700-m intervals on the base map with an arbitrary datum of 53,000 nanoteslas (nT). Several corrections were applied to the magnetometer readings. One was the “recorder drift” amounting to ±0–15 nT. This was determined by periodically comparing annotations on the strip chart of the digital magnetometer display versus the value on the strip chart plot scale. Secondly, the intersections of tie lines with flight lines were used to scale and adjust the diurnal variation values determined from the observatory record. The recorder drift and final diurnal values (amounting to 40–65 nT on the 6 July 6, 1971, flight, and 9–16 nT on the July 7, 1971, flight) were subtracted to yield final magnetic anomaly values and annotated on the base map; the resulting map sheet was hand contoured. The map sheet with the annotated points is still available and has been digitized, rectified, and used to generate the aeromagnetic anomaly grid, so the bias from hand contouring has been avoided. A shaded relief map of the aeromagnetic anomaly is shown in figure 11.

Ground Magnetic Profile

A ground magnetic profile of the vertical field component was carried out by H.R. Blank and M.E. Gettings during August 14–16, 1971. The instrument used was a Scintrex model MF-1-100 fluxgate magnetometer, serial number 809371. The unit and battery pack weighed 3.6 kilograms (kg) and was suspended on a harness in front of the operator. It was necessary for the operator to remove all metal from his person as belt buckles, pens, and so forth, would cause errors in the readings. Readings were taken from the galvanometer on the instrument and recorded to the nearest 5 nT. The ground magnetic profile began in the gneiss of the Uttental Plateau, crossed the marginal border group, all the zones of the layered series, and finished on the Basistoppen sheet (fig. 2), extending approximately 9,020 m in two segments. The location of the profile shown on figures 1 and 2 is approximate because precise station locations have been lost as they were pinpricked on the aerial photos that are no longer available. The interval between readings was approximately 100 m, determined by pacing. A large magnetic offset occurred at 0500 Norway Standard Time, shown by the record (fig. 12) from the Leirvogur Geomagnetic Observatory, Iceland. Because of the highly magnetic rocks of the Skaergaard intrusion, the offset was as much as +2,300 nT, then recovering with time. The drift curve of figure 12 was constructed using the reoccupations of base and some field stations together with the shape of the Leirvogur Observatory record, and the observed data were corrected using this drift curve.

The final corrected profile is shown in plots for the north and for the south segments, (figs. 13 and 14, respectively). In the southern segment, the re-occupations during the large offset can be seen between the 1,000 and 1,500 m distances and show that the drift was satisfactorily compensated for the purposes of this study.

Physical Property Data

Specimen Bulk Density

Bulk density determinations were carried out using the weight-in-air and weight-in-water method. Specimens were soaked in water at least overnight to fill any pores connected to the surface of the specimen. Repeat measurements with the apparatus used yielded a precision of ±0.005 g/cm³ or better for all specimens measured. The sample was weighed three times, together with measurement of the water temperature to obtain the necessary data for determining bulk and grain density,where

w_air

weight of the specimen in dry air,

w_water

weight of the specimen saturated and immersed in water, and

w_airwet

weight of the specimen saturated but surface dry in air.

Then the dry bulk density ρ_bd is given by

${\rho }_{bd}=\frac{w_{air}}{w_{airwet}\mathrm{ }-{\mathrm{ }w}_{water}}{\rho }_{water}$

(1)

where

ρ_water

is the density of water at its temperature during the measurement (Jones and Harris, 1992).

The saturated (wet) bulk density ρ_bw is given by

${\rho }_{bw}=\frac{w_{airwet}}{w_{airwet}\mathrm{ }-\mathrm{ }{w}_{water}}{\rho }_{water}$

(2)

The specimens measured were all non-porous so there was no significant difference in wet and dry bulk density. In the case of volcanic lavas, virtually all vesicles were filled with various zeolites and calcite so that the porosity was near zero. In addition to the 238 measured specimens, specific gravity measurements for 38 specimens reported in Wager and Deer (1939) were included in the dataset in the accompanying data release (Gettings and Parks, 2025). Table 1 gives the distribution of densities for the various lithologies and the units of the Skaergaard intrusion. Figure 15 shows a boxplot of the median densities for units of the Skaergaard and the host rocks. The median density for the gneiss is 2.675 g/cm³ and that of the basaltic lavas is 2.950 g/cm³. Since the northwest half of the intrusion host rock is gneiss (and presumably at depth beneath the basalts to the south) and the southeast half has basaltic lavas as host rock, a Bouguer gravity reduction density of 2.77 g/cm³ was chosen for reduction of the gravity survey data described.

Specimen Magnetization

Eighteen oriented specimens were collected of which fourteen were from the intrusion and four were from basalt dikes, the Basistoppen sheet, and a granophyre dike (table 2). These specimens were core-drilled and measured in a spinner magnetometer (Doell and Cox, 1965) located with the Geophysics Group at Oregon State University. No demagnetization was carried out because the demagnetizing instrument was broken. The magnetic susceptibility of the cores was measured with a susceptibility bridge, also available at Oregon State University. The resulting data were reduced to produce remanent magnetization directions and intensities using the algorithm given in Doell and Cox (1965). Table 2 reports the magnetization properties for these specimens and figures 16 and 17 show the reversed and normal remanent directions on stereographic plots. The normally polarized specimens have a Q (Koenigsberger ratio of remanent magnetic intensity to Earth field induced magnetic intensity, Grant and West, 1965) of 0.09 (granophyre) to 3.91 (basalt dike). With the exception of one specimen (HRB–7A, from the chilled margin in a disrupted zone at the northern intrusion contact), the intrusion specimens are all reversely polarized with an average inclination of −47.7°, declination of 187.2°, and average intensity of magnetization of .003 electromagnetic system of units per cubic centimeter, median 0.00316±0.00582. Median Q is 1.64±0.58 (interquartile range). This corresponds to an average susceptibility of 0.0034 centimeter-gram-second electromagnetic system of units (cgs) with the remanent 1.64 times this in the reversed direction. Subsequent to the work reported here, Schwarz and others (1979) reported their results of samples collected in 1976 using diamond drill sampling and including demagnetization analysis. Their results are an average inclination of −59°, declination of 170°, and they are likely more accurate than those reported here because of uncertainty introduced using oriented hand specimens rather than drilled cores.

Table 2.    

Measured natural remanent magnetization and magnetic susceptibility of oriented specimens, Skaergaard intrusion, East Greenland.

[All intrusion specimens are reversely polarized, except HRB–7A. Remnant magnetization directions and intensities were produced using the algorithm given in Doell and Cox (1965). Abbreviations: ID, sample identification number; I, inclination degrees positive downward from horizontal; D, declination degrees clockwise from north; Jrem, centimeter-gram-second electromagnetic system of units intensity (cgs), units of 10-3 electromagnetic system of units per cubic centimeters; K, cgs magnetic susceptibility, units of 10-6); Q, Koenigsberger ratio; CMGabbro, chilled margin gabbro; cm, centimeter; L a, lower zone unit a; m, meter; L b, lower zone unit b; L c, lower zone unit c; Mississipian Late_<, middle zone; Mt, magnetite; Early_A b, upper zone unit b; Early_A c, upper zone unit c; opx, orthopyroxene; aug, augite; RevAvg, average magnetization direction and intensity of all reversely magnetized samples; NA, not applicable; Is., island; NrmAvg, normal average]

Table 2.    Measured natural remanent magnetization and magnetic susceptibility of oriented specimens, Skaergaard intrusion, East Greenland.

ID	I	D	Jrem	K	Q	Lithology
Skaergaard intrusion
HRB–7	−61.0	261.3	3.52	1,786	3.76	CMGabbro 15 cm from gneiss contact
HRB–8.7	−33.3	226.3	0.87	1,992	0.83	Cross-bedded La
HRB–15	−57.8	209.8	3.16	2,962	2.03	Basalt dike 6 m thick
HRB–15A	−66.9	193.8	0.46	1,135	0.78	L b gabbro
HRB–22	−74.7	181.1	0.46	1,350	0.66	L c gabbro
HRB–34	−52.7	176.6	6.89	13,608	0.97	Mississipian Late_< gabbro Mt rich
HRB–41	−36.3	175.5	8.89	13,589	1.25	Mississipian Late_< gabbro Mt rich
HRB–64	−50.3	166.2	0.06	389	0.30	Basalt dike
HRB–71	−14.6	21.6	0.26	120	4.09	Early_A b gabbro
HRB–75	−59.3	194.7	6.28	3,754	3.19	Early_A c purple band
HRB–84	−40.3	173.0	0.99	226	8.36	Basistoppen sheet, middle zone opx-aug gabbro
HRB–102	−37.4	156.1	3.45	3,246	2.03	CMGabbro Mellemø Island about 1 m from contact
HRB–105	−70.8	198.0	6.56	4,714	2.66	Mt rich layer Nunatak II
RevAvg	−47.7	187.2	3.20	3,739	1.64	NA
HRB–7A1	79.8	176.8	0.85	8,362	0.19	CMGabbro about 1 m from gneiss contact
Samples of non-Skaergaard intrusion rocks
HRB–56A	74.7	54.7	0.06	1,176	0.09	Granophyre vein 0.3–1.0 m thick
HRB–59	63.9	17.5	6.67	3,254	3.91	Basalt dike
HRB–88	68.9	346.2	4.09	6,744	1.16	Dark green px gabbro Basistoppen, midway sta 88 to 89
HRB–101	58.8	12.0	3.19	2,366	2.57	Post-Skaergaard dike 0.7–1.3 m thick, Mellemø Is.
NrmAvg	68.0	16.9	3.00	4,380	1.31	NA

Quantitative Interpretation

Gravity Survey Analysis

Based on geologic mapping and analysis of the differentiation trends (cryptic variation, Wager and Deer, 1939) of the rocks from the exposed zones and the chilled margin of the marginal border group, Wager and Deer (1939) predicted that the intrusion must include what they call a “hidden zone” amounting to 60 percent of the volume of the original intrusion. After further analysis, Wager and Brown (1967) estimated the hidden zone to be 70 percent of their original estimated volume of 500 cubic kilometers (km³) for the entire intrusion, giving a volume estimate for the hidden zone of 350 km³. The cross section given in Wager and Brown (1967) shows the hidden zone as a cylinder of about 7.5 km diameter at the base of the lower zone and narrowing to 5 km diameter at the bottom. Using an average diameter of 6.25 km and volume of 350 km³ this gives a height of 11.4 km for the cylinder. The continental margin flexure, which averages 22° down to the southeast across the intrusion, should make the hidden zone more nearly a vertical cylinder, judging from the cross section given in Wager and Brown (1967). Assuming the contact between the hidden zone and the lower zone to be approximately planar, the contact is at sea level at the northern end of the intrusion and at about 4 km depth at the southern end. The expected gravity anomaly on the axis of such a cylinder, computed from the formulae of Dobrin (1960) with a density contrast of 0.3 g/cm³ would be 30.6 mGal at 0 km depth (northern end of the intrusion), 17.0 mGal at 2 km depth (middle of the intrusion), and 9.6 mGal at 4 km depth (southern end of the intrusion). Calculating a horizontal profile for such a cylinder (Nabighian, 1962) gives an expected width of the anomaly at half amplitude (half width) of approximately 7 km. This estimate does not include the gravity effect of the exposed zones that are below sea level (fig. 3; also fig. 9 in Wager and Brown, 1967). Estimating this effect with half a cylinder of 4 km radius, density contrast of 0.42 g/cm³ (average of lower, middle, and upper zones) and 2-km thickness gives 13.4 mGal at the center of the intrusion; together with the 17.0 mGal from the hidden zone at the middle of the intrusion, the anticipated gravity effect from the Wager model of the intrusion below sea level is about 30 mGal.

A regional gravity anomaly trend surface was computed using orthogonal polynomials (Oldham and Sutherland, 1955) with a least-squares fitting program on the 1 km gravity anomaly grid to account for the large, approximately north-south gradient across the variable density Bouguer anomaly field (fig. 10). The program, originally written by Ron Wahl of the U.S. Geological Survey (USGS), is available as program “surfit” in Phillips (1997). Fits of first, second, third, fifth, seventh, and eighteenth order were examined. The coastal flexure gradient was judged best fit by the fifth order fit (fig. 18) and yielded a residual gravity anomaly map (fig. 19). This residual anomaly map was used for all interpretive modeling. The eighteenth order fit was examined to see if it made a suitable smoothed map but even at this high order too much of the shape and amplitude of the map anomalies were lost.

The residual gravity anomaly map (fig. 19) shows a north-trending, dike-like positive anomaly in the western part of the intrusion outcrop. The anomaly maximum varies from 7 mGal near the northern end to 20 mGal near the southern end of the intrusion, consistent with southward tilt and increasing volume of the intrusion below sea level. The half-width of the anomaly varies from 1.5 km in the north to approximately 3.5 km in the south-central portion to 2.5 km in the south. The maximum in the south occurs beneath the Basistoppen gabbro sheet (McBirney and others, 1989) which may contribute several mGal to the anomaly here and along the extent of the sheet in a northeasterly trend. To the east of the intrusion contact, a broad positive anomaly of approximately 5 mGal is probably due to a combination of the picrites in the Vandfaldsdalen area lava and gabbro including the Macrodyke (M in fig. 2) that cuts through the northern end of Vandfaldsdalen (McBirney and others, 1989). East of the intrusion’s northeast contact, the concentric low of −19 mGal and its extension to the southwest is likely due to the Haengfjeldet formation tuffs (McBirney and others, 1989) that were not accounted for by the variable density Bouguer reduction. The negative anomalies of −4 to −10 mGal to the southwest of the intrusion contact are due to seawater as bathymetric depths were not included in the terrain corrections. The anomaly lows and highs to the west and northwest of the intrusion contact are attributed to lithologic variations in the Precambrian metamorphic complex that contains both amphibolitic (dense) and quartzo-feldspathic (less dense) zones, as well as post-Skaergaard-age basaltic dikes and gabbro bodies (McBirney and others, 1989). Finally, the broad gravity anomaly high of 8 to 12 mGal to the south-southeast of the intrusion contact is due to picrite and dolerite sheets in the basaltic volcanic section (McBirney and others, 1989; Wager and Brown, 1967). The gravity anomaly due to the hidden zone of the intrusion and the below–sea level rocks of the lower, middle, and upper zones is thus assumed to be the dike-like anomaly with a 2–3.5 km half-width and maximum amplitude of 18 to 20 mGal. The possible errors in anomaly values from the various sources is probably about ±2 mGal, so the estimated anomaly maximum amplitude is likely 16 to 22 mGal.

Three-Dimensional Models

The excess mass (Grant and West, 1965) of the intrusion beneath sea level was calculated using a smoothed version of the residual gravity map over a 20×20 km area centered on the intrusion. The smoothed map removed the highs and lows outside the intrusion to reduce error. The integration was carried out using Simpson’s Rule (Abramowitz and Stegun, 1964) and the correction for finite integration area (Grant and West, 1965) was applied. Depending on whether one uses a horizontal cylinder or a sphere approximation to estimate the depth to center of mass, the estimated mass of the intrusion beneath sea level is 13.17×10¹⁵ to 13.58×10¹⁵ g. Because the intrusion is elongate in the north-south direction, choosing an intermediate value closer to the horizontal cylinder gives 13.3×10¹⁵ g as the best estimate.

Several 3D models were computed to further define the sub-sea level shape of the intrusion. First, the program of Cordell and Henderson (1968) was used to model the smoothed residual gravity anomaly grid iteratively with prisms beneath each grid point extending from sea level down to a depth such that the model fits the gravity grid. Two such models were computed, one with a prism density contrast of 0.3 g/cm³, a clear minimum considering the sub-sea level middle and upper zone rocks, and 0.4 g/cm³, a better estimate of the sub-sea level intrusion average density contrast. The 0.3 g/cm³ model is published in McBirney (1996) as a contour map figure. This model has two maximum depths of 3.7 km beneath the central maximum and 3.8 km beneath the southern maximum. The 0.4 g/cm³ model is probably more realistic, at least in the southern half, and has maximum depths of 2.7 and 3.0 km (fig. 20) and is the best estimate of the sub-sea level part of the intrusion using this method.

To make a more realistic model than one with a constant density contrast, the polygonal layer model of Talwani and Ewing (1960) was considered. This model is designed for horizontal layers that would require complicated estimates of the density contrast that would be applied to the Skaergaard intrusion, given its southward tilt. Instead, a more appropriate model is one with dipping polygonal layers of constant density contrast that simulates the average dip and strike of the intrusion zones. To do this, the horizontal components of the gravitational attraction of a polygonal layer were derived by the author using the computer code of Talwani and Ewing (1960) modified to calculate the vector components. Then the gravity effect (vertical component of attraction) of a dipping layer could be obtained by a simple coordinate transformation applied to the gravitational attraction vector for the dipping layer. Examination of geologic map relations (Wager and Brown, 1967) indicates the zones have an average strike of N.55° E. and average dip of 22° SE. The calculated grid of the Cordell and Henderson (1968) model depths (fig. 20) was projected to this plane, and the contours of equal depth in the dipping coordinate system were used to define the polygons for the zones. Density contrasts for the zones were: overlying basalts, upper border group, and Basistoppen sheet, 0.27 g/cm³; upper zone, 0.50 g/cm³; middle zone, 0.47 g/cm³; and lower zone and hidden zone, 0.27 g/cm³. In all, 31 layers were used in the model. The coordinates of the 82 gravity stations in the intrusion and immediate surroundings were also projected so that the model was calculated at the actual altitudes and locations of the stations as opposed to the Cordell and Henderson model that is calculated on a sea-level plane. This is important because most of the stations in the central and eastern part of the intrusion are at high altitudes whereas the stations of the western part of the intrusion are at or near sea level. Figure 21 shows the projected, dipping contour map of the model in figure 20. Figure 22 is a contour map of the final dipping polygonal-layer model gravity effect, including the smoothed residual gravity anomaly from figure 19 for comparison. This model fits the residual gravity anomaly within ±2 mGal in the area within a 1 km buffer about the intrusion contact except in the southeast where dolerite sills in the basalts were not modelled. The mean misfit at the stations is −1.4 mGal and the root-mean-square misfit is 3.67 mGal, due to the negative anomalies in the gneiss to the west and the sediments/tuffs to the northeast. Within a 1 km buffer around the intrusion contact, the mean misfit is −0.65 mGal with a root-mean-square misfit of 3.0 mgGal. The density contrasts could be changed slightly for one or more layers to improve the fit. However, since the fit was well within the ±2 mGal estimated error envelope, no further models were computed.

Because further understanding of the crystallization and differentiation sequence of the intrusion was a primary objective of the research effort, the final calculations were to estimate crudely the relative volumes of the various zones. Subsequent erosion of the intrusion and surrounding host rocks preclude determination of the original extents of the zones, so a tilted-right elliptical cylinder truncated by the horizontal plane of sea level was used as a model for the above–sea level part of the intrusion (fig. 23). This seemed a reasonable model since it is unknown whether the original zones were larger or smaller than a cylindrical model. A best-fitting ellipse was calculated fitting the projection of the intrusion margin dipping at 22° southeast with strike N. 55° E., resulting in a cross-sectional area of 54 km². Appendix 4 details the calculation of the volume of the wedge-shaped portion of the elliptic cylinder. Below sea level, the polygonal layered model discussed above together with the average densities was used to calculate the zonal volumes and excess masses relative to the gneiss density of 2.73 g/cm³. Layer volumes were calculated using the formula for polygonal area (Harris and Stöcker, 1998) for each lamina and then numerically integrated across the thickness of the zone. Figure 23 shows the results of these computations; it is notable that the estimated volume of the hidden zone is only about 1 percent of the intrusion volume.

Three-Dimensional Modeling of Aeromagnetic Data

Interpretation of the aeromagnetic map was focused on the Skaergaard intrusion, although study of the map (fig. 11) shows significant structural boundaries trending approximately east-west across the map area. In part, these boundaries are related to rock units of reverse and normal magnetization, especially the intrusion and the older basalts to the east being reversely polarized against the normally polarized, younger basaltic section to the south, but there are also significant east-west magnetic boundaries to the north in the Watkins Fjord and gneiss terrane (figs. 1 and 2). Across Kangerdlugssuaq Fjord to the west (fig. 1), the Kangerdlugssuaq complex and Kap Edvard Holm intrusive complex 10-15 km southwest of the Skaergaard intrusion (Wager and Brown, 1967) have positive anomalies and so appear to be normally polarized. Within the Skaergaard intrusion the normally polarized Macrodyke (unit M in fig. 2) appears to persist as a positive anomaly that reduces the total amplitude negative anomaly from the Skaergaard intrusion approximately east-west across the intrusion in the subsurface into Kraemers Island (fig. 11). Several faults mapped within and around the intrusion (McBirney and others, 1989) may be correlated with some horizontal gradients on the aeromagnetic map.

The quantitative interpretation was carried out in early 1972, before the 3D Talwani-style polygonal layer software was available. Instead, the author had completed generalizing the right circular cylinder model of Reilly (1969) for an arbitrarily oriented cylinder, rather than one with a vertical axis, and this was used to model a profile across the intrusion. The profile B–B′ was chosen along a north-south line over the highest topography (fig. 11) to obtain the best signal-to-noise ratio at the constant flight height of 1,500 m (5,000 ft). The profile data were taken at 500 m intervals from the aeromagnetic anomaly map and carefully checked against the flight line strip charts to obtain the best estimate of diurnal variation and guard against contouring errors. A linear increasing southward drift of 118 nT was removed from the profile. The resulting magnetic variations (circular symbols in figure 24) were modeled using dipping cylindrical slabs. Radii and centers of circular cylinders approximating the topography and intrusion zones were determined using the topographic map and geologic map of Wager and Brown (1967) and are shown in plan view in figure 25. To approximate the dip of the intrusion zones, the axes of all the cylinders dip at 22° to the southeast, so the circles shown in the plan view are actually somewhat elliptical. In all, 17 cylinders were included in the model, however the layers below the base of the lower zone (base of cylinder 12, fig. 25) only contributed a smooth, long-wavelength offset to the model and do not contribute to the profile model shape at lengths of 10 km and less. Figure 26 shows a cross section of the cylinders model along the profile. Cylinder 16 extended the hidden zone to 5 km depth and cylinder 17 extended the hidden zone of Wager and Brown (1967) down to 20 km depths they proposed. Cylinders 16 and 17 do not contribute significantly to the model which is dominated by the magnetite-rich middle and upper zones. The large contributions of the middle and upper zones make it difficult to define any anomaly due to the hidden zone.

The geomagnetic induction field at the Skaergaard intrusion at the time of survey (July–August 1971) had a declination of 325°, an inclination of 79° downward, and an intensity of 53,500 nT (International Association of Geomagnetism and Aeronomy [IAGA], 2020). The relative magnetite contents of the zones were estimated using the Cross, Iddings, Pirsson, Washington (CIPW) norms (a standard petrologic method to calculate the mineral composition from the bulk chemical analysis, Barth, 1962) from the average zone compositions given in Wager and Brown (1967). These magnetization estimates were combined with the induction vector and the remanent direction from the measurements given above to arrive at an effective magnetization vector of declination 145° and inclination 65° upward. The resulting effective susceptibilities (cgs units) were: upper border group, 0.008; Basistoppen sheet, 0.006; upper zone, 0.010; middle zone, 0.009; lower zone, 0.005; and hidden zone, 0.002. The hidden zone value was calculated from the CIPW norm of the proposed hidden zone composition given in Wager and Brown (1967). Note that these estimates are relative and can be scaled up or down as a group to fit the observed magnetic profile. The upper and middle zones contribute the main part of the model anomaly with lesser contributions from the other zones, as clearly shown on the northern and southern ground magnetic profiles (figs. 13 and 14). The least squares fit of the cylinder model to the observed profile is shown in figure 24 with the plus symbols. Note that the model does not include any contribution from the normally polarized dolerite sill on the south margin of the intrusion (grid north direction 4.5) which would increase the model response between 4.5 and 7.5 grid north direction, and the sharp relative positive peak in the middle profile is inferred to be the superposition of the Macrodyke in the shallow subsurface beneath the profile. The superposition of the Macrodyke may also be present in the northern segment of the ground magnetic profile (fig. 13) between profile distances 2,800 and 4,000 m because magnetite appears as a cumulus phase in large quantities in lower zone C and the middle zone (Wager and Brown, 1967) and would be expected to yield a more negative anomaly. It was concluded that the exposed rocks of the intrusion and their subsurface extent account for essentially the entire aeromagnetic anomaly and further modeling was not attempted.