Introduction

Thawing of ice-rich permafrost and associated massive ice deposits can create severe environmental and engineering problems. The usual solution to these problems are avoidance and specialized construction techniques. Both solutions require a precise knowledge of the permafrost distribution and the position and geometry of massive ice and ice-rich permafrost; this know-ledge depends on the ability to detect massive ice and ice-rich permafrost and to disting-uish between permafrost and unfrozen soils.

The density of permafrost soils in interior Alaska differs markedly from that of segregated ground-ice deposits. For example, permafrost composed of saturated silt with no segregated ice may have a bulk density of about 1.6 M_g m³ compared with a density of about 0.9 M_g m³ (or even less if it contains significant trapped gas) for deposits of ground ice that do not contain soil particles. The moist bulk density of soils sampled in interior Alaska ranges from 1.35 to 1.70 M _g m³ (Furbush and others 1980). Such a relatively large density contrast suggests that massive deposits of ground ice may he detected as an anomaly in gravity measurements. Reference MackayMackay (1962) and Rampton and Malcott (1974) have carried out gravity measurements in the Mackenzie delta region which have shown that the large ice cores of pingos can be readily detected with a sensitive gravimeter.

For use as a method of detection for ground ice in environmental and engineering problems, gravity measurements must be capable of detecting much

smaller ice masses than are found in most pingos. The procedures and analysis of Rampton and Walcott (1974) suggest that relatively small ice masses may be detectable provided that the spacing between measurements is reduced and a very sensitive gravimeter is used,

This paper reports gravity measurements made over deposits of ground ice that were relatively small compared to pingos and for which there was no appar-ent surface expression in the form of topographical or vegetative changes. The measurements were carried out along a profile over undisturbed terrain at a site where a road cut was later made. Ground truth was obtained by reconstructing a profile of subsurface conditions at the site from pre-construction borings and from information gathered at the time that the cut was made (Osterkamp and others 1980). After completion of the roadbed, a second set of measurements was made along one edge. A third set of measurements was made across an artificially-constructed ice mass of known dimensions (Gruol and others 1981).

Experimental Procedure

The main study site was the road cut at Engineer Creek, located about 9 km north-north-east of Fairbanks, Alaska, along a realignment of the Steese Highway. The site lay on a gently-sloping, north-facing hillside at elevations between Z71 and 282 m. The first measurements were made between survey stations 462+00 and 466+00 (State of Alaska, Department of Transportation and Public Facilities designation), which are 122 m apart. In the study area, there was ~5 to >12 m of silty soil underlain by bedrock. Small black spruce trees, generally <10 m in height, covered the entire area. Further details on this site have been reported by Osterkamp and Ourick (1980).

The gravity stations were surveyed at intervals of 3.05 m along a line displaced 24 m left (facing north-west) of the centerline. Plywood boards 0.3 m in diameter by 16 mm thick were placed on the mineral soil and fixed in place by spikes of length 0.2 m driven into the soil. Elevations of these boards were measured to a precision of ± 6.1 mm. The gravity measurements were marie with the gravimeter resting on the boards.

The instrument used was a Lacoste and Romberg gravimeter (Model G, S/N 248). Accuracy is stated by the manufacturer to be ± 0.01 rnGal (1 Gal (Galileo) = 1 cm s ² = 10 ² m s ²). Repeated measurements over a test station indicated the actual sensitivity was closer to ± 0.004 mGal, although such sensitivity is rarely obtainable under field conditions.

Gravity profiles were obtained on 2 and 5 August 1976 proceeding from station 462+00 to 466+00. Each set of measurements was completed in ~ 4 to 4.5 h. During the course of the profiling, three measure-ments were made at station 462+00 which were used to correct for drift and earth tides. The maximum possible difference from linearity for the earth tides over a period of 2 h was less than 0.004 mGal during the two days so that tide and drift errors were considered linear in the subsequent analysis.

After completion of the roadway, a second gravity profile was obtained on 29 October 1979 between stations 454+00 and 463+00 on one edge of the roadbed constructed in the road cut. This profile was made parallel to the first profile hut displaced by about 36 m.

A gravity survey was also made over an artificially-constructed ice mass (Gruol and others 1981). The ice mass of known dimensions and burial depth (34 × 0.69 x. 3.2 m buried at a depth of 1.2 m) was constructed to evaluate various devices for detecting ground ice. The results of this third gravity survey are also reported in the next section.

Results and Analysis

It is usual, in reducing gravity data, to find the Bouguer gravity corresponding to the observations. Accordingly, an arbitrary datum plane was chosen at the level of station 466+00 (271 m in elevation) and the data and the initial analysis discussed in this section are referred to this plane.

Ho drift/earth-tide corrections were applied to the set of data taken on 2 August 1976, because the drift appeared to be within the nominal precision of the gravimeter. However, the second data set was cor-rected for drift tides (by linear interpolation) with reference to the three readings taken at 462+00. For this second set, the maximum meter-reading drift was 0.028 mGal corresponding to an actual calibrated dif-ference of 0.031 mGal. The average values of the two data sets after correcting the second set of measure-ements for drift and tides are shown in the second col-umn of Table I together with the elevations relative to station 466+00 and the corrections for latitude and free air. Also shown are the readings corrected for drift/tides, latitude and free air in the final column. Since only relative values of gravity were of interest to determine whether it might be possible to detect ice masses along the profile, no attempt was made to tie the survey to absolute gravity values.

To calculate the correction (Ag|_) to the gravity due to latitudinal change in Table I, the relation

was used, where r is the Earth’s radius at the equator, △S is the arc length along the profile, ф is the latitude and θ is the great circle bearing of our study line with respect to true north (Grant and West 1965). For the site at Engineer Creek, ф is 64°55’ and θ ≈ 43°, giving △g_L = 46 × 10⁴ △S mGal where AS is measured in meters proceeding along the study line from north to south,

The free-air correction △g_F- is given by the relation

where Ar is the elevation from the reference level measured in meters.

Note that although the uncorrected values tended to increase with a decrease in elevation along the profile, the values corrected for latitude and free air in the last column of Table I, on the other hand, decrease with decreasing elevation.

The Bouguer correction, which accounts for the underlying material, is given by

where G is the universal gravitational constant, pbis the Bouguer density and h the elevation above some reference plane which is here taken as the horizontal plane through station 466+00. With h measured in meters and pb in Mg m-_g gb = 0.0419 mGal.

The bedrock underlying the Fairbanks region has a density typical of continental crustal rocks. When p₀ is 2.67 M_gg m³, △gh at 462+00 relative to 466+00 is approximately 1.26 mGal. It is clear from this and the other elevation values in Table I that the decreasing trend would not be removed by the Bouguer effect. Figure 1 shows a plot of the residual gravity values along the profile after application of the Bouguer correction; also included is the elevation as a function of station along the profile.

Fig. 1. Elevation h and gravity g (corrected for P = 2.67 M_g m³) latitude, free air and Bouguer terms, as a function of station positions at the undisturbed site at Engineer Creek.

Inspection of the Fairbanks quadrangle map (0-2 NE) and maps of adjacent areas indicates that terrain corrections are negligible in comparison with the gravity anomalies obtained near the road cut at Engineer Creek. Furthermore, the major effect on the gravity profile from the terrain would be to put a linear trend in the data. However, for completeness, terrain corrections at the endpoints of the profile were estimated from Hammer rings {Reference HammerHammer 1939) up to a radius of 4 469 m (rings 3 through J) using the tables compiled by Reference BibleBible (1962). For station 466+00 the terrain correction was approximately Û.25 mGal and for 462+00_approximately 0.26 mGal with pb = 2.67 Mg m 3. As one would expect, the difference (0.26 - 0.25 = 0.01 mGal) arose mainly from the terrain in the innermost rings (within 170 m).

It is interesting to note that, even if all corrections based on surficial features are made, there remains a significant linear decreasing trend with decreasing elevation. In the absence of know-ledge of the stratigraphy of the area, a trend of this type, appearing in a large-scale gravity survey, is usually ascribed to large-scale regional differences in the density of the material underlying the area. Here, however, the scale size of the survey is considerably smaller and, further, the depth to bedrock (schist) was known along the profile at several points from bore holes made by the Department of Transportation and Public Facilities, State of Alaska. The depth to bedrock from the undisturbed surface varied from 4.9 m at station 458+00, 9 m at station 462+00 and >12 m at station 466+00. The bedrock was overlain by permafrost consisting of silt and ice.

It is clear that the sloping (relative to the soil surface) contact between silt and bedrock should affect the gravity profile by tending to cause a lower reading at the lower point (466+00) than at the upper point, i.e., the apparent density below 466+00 is less than that below 462+00.

Since the sloping bedrock is well below the surface and extends considerably beyond the limits of the survey profile, its effect can be considered part of the linear trend in the data (Fig.l). In addition, because the primary emphasis of this report is the use of gravity measurements made for the purpose of detecting deposits of massive ice, linear trends over the dimensions of the profile are not of interest. Consequently, further refinement of the data in its original form in terms of terrain and Bouguer corrections using varying densities has not been undertaken. Instead, linear fits were made to the elevation data and to the gravity data corrected for drift/tides, latitude and free-air terms as tabulated in Table I. Deviations from the least squares fit were then calculated. Eliminating the linear trend should yield gravity deviations from the least squares fit which may be mainly ascribed to local variations in the density of the topmost layers of material which consists mainly of silt and ice.

Some of the variation, however, may arise from undulations of the terrain in the vicinity of the profile line. Having eliminated the free-air gravity before obtaining the deviations in gravity from the least squares fit, only the Bouguer and terrain gravities need be considered. Since only relative

gravity differences are of interest, we may remove the effects of the Bouguer gravity using the deviations from the least squares fit in elevation (Fig.2.: top panel). It should be noted that normally one uses the same datum level for both the free-air and Bouguer gravities, but here we take a new datum level, i.e. the linear least squares fit in elevation, as the datum level to correct for the Bouguer gravity. The results of this procedure are shown in the middle panel of Figure 2. Here we have used a Bouguer density of 1.6 Hg m”3 for the unfrozen, moist, surface soil layer, chosen as representative on the basis of a few bore-hole samplings available from the site and also on a number of soil samplings made in interior Alaska (Furbush and others 1980). (Note that, strictly speaking, one should consider only the component of the vertical field perpendicular to the least squares fit for elevation; however, since the grade is less than 10%, this component differs from the vertical by only about 5 parts in 1 000 so that little error should be introduced by this procedure).

Fig. 2. Deviations from linear least squares fit of elevation (top panel) and gravity (middle panel) first corrected for latitude and free-air and then a Bouguer (with pb = 1.6 M_g m³) term using the elevations deviations. The bottom panel shows the ice masses and other ice-rich areas revealed during road-cutting at Engineer Creek.

The terrain correction from the new least-squares fit, elevation, datum plane for each point must be fairly small and the relative differences for all points even smaller. Along the profile, the deviation in maximum elevation is about 0.7 m (Fig.2.) and while the deviations transverse to the profile are not known exactly, due to lack of a contour map of finer scale than 3.05 m, the gentle slope suggests even fine-scale terrain corrections are not needed. Indeed if the elevation profile were to extend infinitely in the transverse direction, the maximum terrain correctiort would be less than 0.003 mGal (Telford and others 1976: table 2.1 using p_b = 1.6 M_g m³). Thus, the residual gravity profile in Figure 2 must arise mainly from density contrasts in the underlying material within a few meters of the surface but above the bedrock.

Variations in bedrock density may be eliminated as the major source of the residual gravity profile by considering a sphere of radius a = 5 m with center buried at R = 15 m (the bedrock is at least 10 m below the stations). At the surface, directly above it, the sphere will produce a vertical gravity field of g = 1.33 a₃ п G p/R² and the points on the surface where the vertical field drops off to half this amplitude are located at (1)where f, the distance from the center of the sphere to the point on the surface, is defined by r³ = 2R₃. For R = 15 m, the half-width distance is - 23 m which is about the same as the observed half-width for distances between minima in Figure 2. The amplitude of the field of the sphere is 0.0155 p mGal, so that to obtain an amplitude between minima of 0.04 mGal, one needs a density contrast of greater than 2 Mg nr³.Any other mass shapes or shallower burial depth still within the bedrock will require an even greater and more unrealistic density contrast.

The occurrence of massive ice and ice-rich permafrost was documented with pre-construction borings and observations made during road cutting by Osterkamp and others (1980). A scheme of the locations and shape of the ice observed down to the base of the cut is shown in the bottom panel on Figure 2. It is clear that the negative deviations in the gravity profile closely match the locations of the deposits of massive ice.

Under certain simplifying assumptions, a rough estimate of the density contrast along the survey line can be made from the residual gravity profile shown in Figure 2. Since the vertical gravity component only along a single direction was available, it was assumed that the calculated residual deviations were caused by blocks of material with rectangular cross-section with different but uniform densities. Each block was taken to be 3.05 m in width along the profile, 6.1 m from top to bottom and infinite in length transverse to the plane containing the station line and the vertical direction. Figure 3 shows schematically the configuration of the blocks along the profile and the residual (relative to a mean value) deviations that such a collection of blocks might cause. The depth from the surface to the base of the block was chosen as 6.1 m instead of a shallower depth for two reasons: first, because the density contrast would have to be unreasonably larqe to cause the observed residuals, and second, because it is the excess segregated ice in the first 6.1 m below the road bed which may cause problems after road construction. The distance between station points on the survey was 3.05 m so that the block width was taken to be the same. The blocks were taken as infinite in the direction transverse to the profile as a first approximation, since the massive ice deposits appear to be quite long in the transverse direction at this site (Osterkamp and others 1980).

Fig. 3. Scheme of configuration of blocks along profile line used to estimate density contrast in first 6.1 m of depth.

Figure 4 shows the drop-off of the vertical component of the gravity field at the surface from a single block of infinite length, 3.05 m width and depth to the base of 6.1 m and with p = 1 M ₃ m- ³ in a direction perpendicular to the infinite axis. The vertical field at a point y along the surface was calculated using the equation

where G is the universal gravitational constant, p is the density, a is the half-width of the block (1.5 m) and h is the depth to the bottom (6.1 m).

For an infinite slab of uniform density (1 M ₃ m- ³) and thickness (6.1 m), the field at the surface is 0.256 mGal (g = 2TrG_Ph). An infinite slab can be thought of as consisting of an infinite number of blocks with each of them having a vertical gravitational field component of the form of Equation (4). Figure 4 indicates that more than 90% of the field at the central point on the surface can be accounted for by the mass within 21 in to both sides of the central point. In order to make more tractable the problem of inverting the observed field variation along the profile so that the density contrast is obtained, it is assumed here that only the first three blocks to either side of the central point (the point from where data were taken) plus the central block contributed to the observed field. For uniform density, the net contribution from the seven blocks (three to each side plus the central block) accounts for about 83% of the field produced by an infinite slab. (Note that, in fact, the blocks will have different densities in the real situation, but that the density contrast from the mean for a soil-ice mixture should not be very large.) For a given measured field, if it is assumed that only the seven closest blocks have variable density while the blocks outside the centrally-located blocks have a density representative of the mean for all blocks including the central ones, the field at the jth block will be approximately given by

or

where a, b, c and d are geometric factors numerically equal to the values given at positions 0, 3.05, 6.1 and 9.15 m on Figure 4, namely 0.0974, 0.0348, 0.0145 and 0.0076, respectively. If we set g_j - 2_пGph = Δg_j and say ρi – ρ = αi, where p is the mean density, then the above equation becomes

Fig. 4.

Here Ag_-1 may be represented as the deviation of the gravity value from the mean, i.e. the residual gravity values shown in the middle panel of Figure 2, and a-j as the density contrast relative to the mean value of density.

Since no data are available outside stations 462+00 to 466+00, it is assumed that blocks to the left of station 462+00 and to the right of station 466+00 (see Fig.3.) have a density equal to the mean, namely p, so that ak 0 for k = -2, -1, 0, 42, 43, 44. Thus, if the △gj of Equation (5) is taken to be the 41 residual values shown in Figure 2, there will be 41 equations in 41 unknowns.

The set of 41 equations has been solved by least squares using the Gauss-Seidel iterative method. Figure 5 shows the resulting density contrasts plotted as a function of station number. Comparison of the plat of density contrast with that of gravity in Figure 2 shows that the two are quite similar.

Un estimate of the excess ice in the blocks may be made from the calculated density contrasts assuming reasonable values for the densities of the ice, soil and mean density for the layer which is 6.1 m thick. If we take the soil density to be 1.6 Mg m~3, the ice density to be 0.92 Mg m⁻³ and the mean density to be 1.45, the fraction of excess ice by weight is

where s refers to the soil, I to the ice, V to volume and pj, the density of the jth block. Note that the soil may contain ice and water interstitial to the particle grains, but these components are not considered as excess moisture in this calculation. With the mean density taken as p, pj = aj + p. For aj = -0.4, -0.3, -0.2, -0.1, 0 and +0.1, the excess ice by weight according to Equation (6) is 243%, 113%, 61%, 33%, 16% and 5%, respectively, corresponding to 80%, 60%, 51%, 37%, 22% and 7%, respectively, by volume. Clearly the 80% value for the volumetric fraction of excess ice for the block at 462+30 (see Figs.2 and 5) is too high; however, the other values appear reasonable for the excess ice content as can be inferred from Figure 2 (see also Osterkamp and others 1980: 279). No direct integrated averaged values for moisture content are available.

Fig. 5. density contrast (M_g m³) versus distance along the profile line obtained by inverting gravity data.

A second gravity survey was made along the edge of the roadway after paving. Unfortunately, the profile line which extended from stations 454+00 to 463+00 was surveyed only once and some of the values are suspect. However, because some ground truth was available and also because the survey was conducted over a fairly even surface with a gently sloping grade, the data were partially analyzed and are discussed here.

Figure 6 shows the gravity values after drift, tides, latitudinal and free-air corrections were made. The survey was made over a period of 6 h on 19 October 1979, with station 454+00 used as the reference point to monitor the drift and tidal effects approximately once every hour. The excavation to grade required cutting into bedrock between stations 454+00 and 460+55 and subexcavating through silt and ice .between 460+55 and 463+00. As indicated in the elevation scheme shown below the gravity profile, the sub-excavation was backfilled with thaw-stable materials.

Fig. 6. Elevation h and gravity g (corrected for latitudinal and free-air terms) versus distance along one edge of completed road at Engineer Creek. Approximate shape and location of backfill area are also shown.

A linear least squares fit was made to the data between stations 455+00 and 459+00 where there was no obvious large deviations which might be erroneous, and also to avoid edge effects from the backfill material whose bulk density is significantly lower than the bedrock material (schist). The least squares fit is shown as a solid line in Figure 6 extrapolated to the ends of the profile line. It is interesting to note that the backfill shows up in the profile as a region of mass deficit as might be inferred from its density which is about 2.2 M_g m³ in contrast to the bedrock which should be approximately 2.67 M_g m³, typical of unweathered schist.

The density of the backfill can be estimated from the observed gravitational deficit relative to the value extrapolated from the linear least squares fit which is representative of the density of the schist below the region between stations 455+00 and 459+00. The vertical field due to the backfilled material can be approximated by an infinite slab for points distant from its edges, with the same approximation holding for a slab of the same thickness between stations 455+00 and 459+00. The difference is thus 2irGh(pi - p2> = △g. If △g is taken as the difference (approximately 0.175 mGal) between the extrapolated line and the last value along the profile (that at station 463+00), h taken as the thickness of the backfill (~ 5.5 m) and pi taken as the density of the schist (= 2.67 M_g m₃), we obtain the density of the backfill material p2 1.9 Mg m~3. This is in fair agreement with the actual value of about 2.2 M_g m₃ and is reasonable in the light of uncertainties such as the possibility that the subexcavation may have not been made to bedrock between stations 460+55 and 462+00 as originally designed.

The residual deviations relative to the least squares fit along the profile between stations 455+00 and 459+00 are shown at the top of Figure 7. Unlike the residual profile in Figure 2, no attempt has been made to make a Bouguer correction due to the deviations from a least squares fit to the elevation shown at the bottom of Figure 7, because the deviations are extremely small, less than 0.05 m. Such elevation deviations can at most produce about a 0.005 mGal correction to the residual deviations in the vertical gravity component.

Fig. 7. Elevation and gravity deviations from a linear least squares fit of data between stations 455 and 459 along the edge of the road at Engineer Creek.

On comparing the variations in the residual gravity in Figures 2 and 7, it is seen that the amplitudes are of the same order. Aside from possible errors in measurement, the variations seen in the roadside survey must be attributed to density variations in the underlying material. This is possible since schist is known to vary considerably in density. Also it is likely that there is a varying thickness of the aggregate base and sub-base and perhaps even residual pockets of soil overburden which are reflected in the residual gravity profile.

A third gravity survey was carried out in a direction perpendicular to the long axis of an artificially-constructed ice mass with approximate dimensions of 34 × 0.69 × 3.2 m buried to a depth of 1.2 m to the top of the mass (Gruol and others 1981). The ice mass is located in a flat area of stunted spruce where bore-hole, resistivity and other geophysical investigations suggest that there is no naturally-occurring massive ground ice. There is a somewhat ice-rich layer at a depth of about 1 m, but all other data indicate relative lateral and vertical (to a depth of ~ 18 m) homogeneity of the silt.

It is clear that as one approaches the dike of ice, its gravitational field should approximate that of an infinite prism with rectangular cross-section. An equation similar to Equation (4) can be obtained for the vertical gravity component valid along a profile perpendicular to the strike of an infinite prism buried at an arbitrary depth. With a density contrast of 0.68 M_g m³ (p = 0.92 M_g m³ for ice and p = 1.6 M_g m³ for the surrounding silt), the maximum field component deviation along the profile relative to a uniform background field is calculated to be about 0.008 mGal.

If the dike extends to the surface with the same density, the field deficit is 0.022 mGal. Because the average density of the backfill over the ice mass is expected to be somewhat less than the surrounding undisturbed terrain, the maximum vertical field deviations should lie between 0.008 and 0.022 mGal. Figure 8 shows the results of the gravimeter profile across the ice mass near its center together with relative station elevations from a reference plane. Most probably, due to difficulties in stabilizing the supports for the station platforms, the readings clearly show significant variation, even though they have been corrected for latitude, free air and Bouguer terms. However, there appears to be a reduction in gravity directly over the location of the ice mass relative to the four values on either side. This reduction is about 0.02 mGal from the mean of the eight values.

Fig. 8. Elevation and gravity profile data over the long axis of the artificial ice mass.

Conclusions

The usefulness of gravity measurements for oil, gas and mineral exploration is well-established. Aside from the relative scale, the detection of ground ice by means of gravitational measurements is similar. The major differences that exist between exploration of resources and detection of ground ice by gravity arise from the localized nature of the deposits of ground ice, smaller size of the ice masses of interest to the construction industry and the correspondingly greater requirement on sensitivity.

The present investigation indicates that it is possible to detect isolated, massive ground ice deposits by means of commercially-available gravimeters. Such detection is predicated on the size and depth of burial of the ice masses and the density contrast between ice and the surrounding soil. As an example of an ice-mass configuration that is on the threshold of detectability, consider a disk with a diameter of 2 m with circular faces parallel to the horizontal plane, buried to a depth of 1 m to its top face and with a thickness of 1 m and density contrast of 0.68 M_g m³. The contrast in gravity can be calculated to be 0.012 mGal, which is close to the detection limit of commercially-avail able gravimeters. Also it has been shown that an artificially-constructed ice mass of dimensions 34 × 0.69 × 3.2 m buried at a depth of 1.2 m to its top face cannot be detected with certainty by a commercial Lacoste-Ramberg gravimeter.

The use of a commercially-available gravimeter to detect smaller ground ice masses, small ice lenses and ice-rich soils does not appear to be possible, unless they are closely packed together so that they form a distinct gravitational signature. Further, sensitivities as low as 0.01 mGal or less are not attainable without an extremely stable gravimeter platform requiring tedious leveling adjustments.

Without commercially-available field instruments with sensitivities well below 0.005 mGal and with quicker response time from set-up to reading, the gravity technique is not suitable for rapid reconnaissance of permafrost terrain. However, the technique appears to be feasible for detailed work such as a pre-construction investigation of building sites where massive ground ice may occur in an overburden which is otherwise uniform or for investigating selected sites along airfields, roadways and other linear constructions. Estimates of the excess ice content from the gravity data may then be made using analyses similar to the one outlined in this

report or using those widely available to the geophysical exploration industry. At present, the gravity method and possibly the impulse radar are the only non-contacting means of obtaining information on the content of excess ice in permafrost.