Subglacial seismic reflection strategies when source amplitude and medium attenuation are poorly known

Abstract

Seismic reflection techniques are a powerful way to probe physical properties of subglacial strata. Inversion of seismic data for physical properties may be hampered, however, by lack of knowledge of the source amplitude as well as lack of knowledge of the compressional and shear attenuation in the ice. New methods are described to measure the source signature that require no a priori knowledge of the ice attenuation profile. Another new method is described to obtain the angular dependence of the subglacial bed reflection coefficient that is relatively insensitive to knowledge of the ice attenuation. Finally, a correction is provided to a long-standing error in the literature regarding measurement of the bed normal incidence reflection coefficient.

1. Introduction

The mechanisms for rapid basal motion of glaciers and ice streams are poorly understood, but vitally important to better quantify the mass balance of Antarctica and Greenland. Internal deformation of the ice is strongly dependent on the crystal orientation fabric within the ice, with strongly oriented ice relatively softer in shear. The presence of basal water in films or in massive units is important both for dynamic glaciology and for locating and characterizing subglacial lakes that might contain extremophilic life.

Seismic wave speeds in ice and in the subglacial sediments and rocks depend on the bulk moduli, the strength of anisotropy and on the densities of the materials. Shear wave velocities are extremely sensitive to the shear modulus and to fluid content, which is related to porosity. Seismic reflection profiling can be used to map gross subglacial structure and to determine properties of the bed, including till porosities (Blankenship and others, 1987). Seismic reflection amplitudes are sensitive to boundaries where acoustic impedance (defined as the product of density and seismic velocity) changes, and to the roughness of the boundary. These parameters, which can be efficiently determined by seismic profiling – rock and sediment type and roughness at the subglacial boundary; fluid presence, fluid layer thickness and porosity of the sediments; thickness of the subglacial sedimentary layer – are all crucial parameters in understanding basal sliding, subglacial till deformation, internal deformation and subglacial hydrology. Indeed, aside from drilling to the bed, which is extremely expensive, seismic imaging is the only way to determine these properties.

The polar glaciological community has successfully exploited travel-time variability at non-normal incidence (Bentley and Chang 1971; Kapitsa and others, 1996) as well as normal-incidence reflection amplitude of the basal reflector to gather considerable information about the properties of that interface in various settings (Blankenship and others, 1987; Smith, 1997; Anandakrishnan and others, 1998). There have also been attempts to exploit the non-normal-incidence reflection amplitude to determine subglacial properties (Nolan and Echelmeyer 1999; Anandakrishnan 2003).

However, several challenges exist that tend to make quantitative estimates of the properties at the subglacial boundary difficult. One challenge is that the compressional and shear attenuation in the ice may be highly variable and poorly constrained. Another challenge is that the source amplitude is often unknown. Both of these limit the ability to quantify the bed reflection coefficient, or in other words lead to bed reflection coefficients with large uncertainties. The large uncertainties have significant implications for understanding the glacial environment and glacial dynamics since uncertainties can preclude distinguishing, say, unconsolidated sediment from crystalline rock or preclude distinguishing deforming from non-deforming sediments.

The objective of this paper is to address both causes of uncertainties and provide a method for obtaining an estimate of the bed reflection coefficient in situations where lack of knowledge of the ice attenuation and/or the source signature would otherwise lead to significant uncertainties. Though the techniques presented here are general, calculations are provided for glacial environments to demonstrate the methods, inasmuch as seismic reflection data analysis in glacial environments is particularly prone to uncertainties in attenuation, and source amplitudes (typically explosives) are also difficult to estimate.

After a brief review of a standard technique for obtaining the bed reflection coefficient (section 2), methods are proposed in section 3 for obtaining source amplitude independent of knowledge of the attenuation. In section 4, another method is proposed for obtaining the angular dependence of the reflection coefficient or the relative amplitude vs angle (rAVA). This method requires an estimate of the medium attenuation, but is substantially less sensitive to attenuation uncertainties than other techniques.

Consider a seismic reflection experiment (see Fig. 1) with source and receiver placed on the surface with the object of estimating the properties of the horizon labeled ‘bed’. The material properties between the source–receiver and the bed can be variable in the z dimension. With sufficient signal-to-noise ratio, both the primary reflection and the first multiple are recorded at receiver (∇) with amplitudes A ₁ and A ₂, respectively, given by

(1)

where A ₀ is the source amplitude relative to a reference range d ₀ = 1 m, H is the thickness of the intervening medium, d is the distance along the path, the γ_i include spreading losses and the effect of the interface on the received amplitude (e.g. for a receiver on a free surface, at normal incidence the received amplitude is double that of a receiver far from the boundary), r is source–receiver offset, R is the reflection coefficient and α is the attenuation. Note that the angles θ_{i b} are measured between the interface and the ray (‘grazing angle’) rather than the normal and the ray (‘incidence angle’).

Fig. 1. Seismic reflection experiment geometry. In the text, the origin of the coordinate system is placed at the location of the source (⋆).

The assumptions are as follows, many of which have been made solely for clarity in presentation. Thickness H is assumed constant over the area illuminated at the bed, and bed roughness is assumed small, which can be formally stated as 2(kσsinθ)² ≪ 1, where k is the wavenumber and σ is the root-mean-square (rms) bed roughness. It is also assumed that the reflection coefficient at the upper boundary is −1 (though including the exact reflection coefficient at the upper interface is straightforward) and that layering in the medium has a negligible impact on the received amplitudes (adding layers does not change the conclusions, only the mathematical complexity). Also, the distance along the ray path has been approximated here by straight-line rays which may be adequate over some angular range, but could easily be computed if the velocity profile is estimated.

It is clear that if observations for only a single bounce (A ₁) are available, then there are three unknowns near normal incidence: source amplitude, A ₀, the bed reflection coefficient, R, and the ice attenuation, α (it is assumed that the spreading loss and thickness, H, are known or can be estimated reasonably well). By combining the first and second bounce amplitudes, a common practice is to reduce the number of unknowns to two. This approach is briefly described in section 2. In sections 3 and 4, a (fortuitous) opportunity is exploited to eliminate two unknowns in the above equations.

2. Bed Reflection as a Function of R, α

A number of researchers (e.g. Smith, 2007) have used the ratio of the amplitudes to make inferences about properties of the subglacial bed. Taking the ratio of the amplitudes of the first and second bed reflections, the unknown source amplitude, A ₀, is eliminated by recognizing that near normal incidence the reflection coefficient is nearly constant with angle:

(2)

In the limit of grazing angle approaching π/2, γ ₁/γ ₂ = 2 since the spreading loss due to distance is simply a factor of 2 and the ratio of other factors (ray tube cross-section and interface effects) is unity. Thus, the normal incidence reflection coefficient becomes

(3)

and, expressing this in terms of energy,

(4)

Note that a number of papers (e.g. Röthlisberger, 1972; Smith, 1997, 2007; Vaughan and others, 2003) in the glaciology literature have a similar expression, but differing by a factor of 1/2 in the exponent (R ² ∝ e^2αH ). The factor of 1/2 originates from Röthlisberger (see his equation (73)) who defined an ‘attenuation constant’ relative to energy rather than amplitude and hence the factor 1/2. Attenuation, however, is fundamentally and very widely defined relative to amplitude, not energy (e.g. Sheriff and Geldart, 1995). The confusion from different quantities with identical units leads the latter papers to erroneously use the amplitude attenuation coefficient in an expression derived from energy considerations. This led to underestimates of the true reflection coefficient and hence estimates of impedance of the subglacial bed.

The magnitude of the errors depends on the ice thickness and attenuation of the particular case. As a specific example, consider a case with H ∼ 2200 m and α = 0.21 × 10⁻³ m⁻¹ (akin to Rutford Ice Stream, West Antarctica; Smith, 1997). In this case, the reflection coefficient was underestimated by a factor of 0.63 relative to the true reflection coefficient given by Equation (4). Thus, the largest reflection coefficient in that study should have been |R| ∼ 0.35 instead of the reported |R| ∼ 0.22. The errors in |R| lead to errors in estimating the acoustic impedance of the bed. The error in bed impedance is a function of |R|, with errors increasing with increasing |R|: for example, for |R| = 0.1, the bed impedance was underestimated by a factor of 0.93; for |R| = 0.3, by a factor of 0.79; for |R| = 0.5, by a factor of 0.64.

Of particular relevance to this paper is the fact that the uncertainties (in |R| and bed impedance) are much larger than were predicted. For example, a normal incidence reflection coefficient of 0.2 would actually have lower/upper bound uncertainties of 0.11 and 0.6 (using the same attenuation uncertainties as used by Smith). The resulting acoustic impedance, 5.2 × 10⁶ kg m⁻² s⁻¹ then has lower/upper bound uncertainties of 4.3 × 10⁶ and 13.9 × 10⁶ kg m⁻² s⁻¹, resulting in a wide range of possible subglacial environments (the true uncertainties in this case are nearly four times larger than reported). The large uncertainties mean, for example, that it can be much more difficult to determine with reasonable confidence whether the bed is deforming or non-deforming.

The main point of this section is that a significant limitation of this method (Equation (4)) is the large uncertainty associated with the unknown (and highly variable) attenuation. Attenuation is a strong function of temperature, which is generally poorly known in glaciers. These large uncertainties lead to difficulty in determining the bed properties and were the motivation for the new techniques described in the following sections.

3. Source Calibration

3.1. Multiple bounce method

Generally, with two equations and three unknowns, one of the unknowns can always be eliminated. However, in the case of Equation (1), one can exploit the exponential factor to eliminate two of the unknowns by taking the ratio

(5)

and, in the limit of θ approaching π/2, d ₂ = 2d ₁ = 4H, a simple expression for the source amplitude is obtained,

(6)

It is important to note that this expression requires no knowledge of either the attenuation or the bed reflection coefficient. Thus, it is an extremely useful way to measure source characteristics and quantify shot-to-shot variability. In an isovelocity medium bounded by a free surface, 1/2γ ₁ = H/2d ₀.

At angles close to normal incidence, the source amplitude can be obtained by

(7)

The accuracy of Equation (7) away from normal incidence depends upon the degree to which the reflection coefficient differs between angles θ _1b and θ _2b. In many cases it is expected that Equation (7) should be valid in the range within approximately 10° of normal incidence.

In addition to being a convenient expression for estimating the source amplitude, Equation (7) also has the virtue that it is relatively insensitive to the uncertainties in the attenuation. Note that the quantity 2d ₁ − d ₂ in the exponent is very small and is identically zero at normal incidence. At angles near normal incidence, the exponential is well approximated (via Taylor series expansion) by exp(3αr² /8H) and in many instances will be nearly unity.

This method provides source calibration only for source angles close to normal incidence, which could be problematic if the source beam pattern were strongly asymmetric. However, in polar environments, the source is often initiated in firn which has a strong positive seismic velocity gradient, leading to near-vertical ray paths for bed-reflected arrivals. Thus, even in an experiment with a source exhibiting a strong angular dependence, the above technique may suffice.

Finally it is noted that if the source has a shear wave component, the source shear amplitude could also be obtained using Equations (6) and (7) with the amplitudes of the first and second shear multiple.

3.2. Direct path method

Here we obtain an expression for the source amplitude using the same strategy as in section 3.1 but using direct path measurements instead of bed reflections. Consider the geometry in Figure 2, where we assume paths that do not interact with the bed boundary. The amplitudes for the direct path are

(8)

where s is path length. Now the uncertainties of the unknown attenuation can be minimized by taking a pair of arrivals where the ratio of the ray-path distances s ₂/s ₁ = 2. If it can be assumed that the depth-averaged attenuation is approximately the same for the two paths (which would be true, for example, if both paths sampled all or most of the firn layer) then the source amplitude is simply

(9)

So whilst knowledge of the velocity profile is required, there is no requirement for knowledge of the attenuation to obtain source amplitude.

Fig. 2. Seismic reflection experiment geometry and ray paths for direct path experiment.

4. Bed Reflection Coefficient Vs Angle

This section is concerned with estimating the angular dependence of the reflection coefficient (AVA) which generally contains useful information about the bed properties.

4.1. Direct estimation of the reflection coefficient, AVA

Given a dataset with one and two bounces from the bed near normal incidence, one could use Equation (6) to estimate the source amplitude (and its shot-to-shot variability) and then rearrange the A ₁ from Equation (1) to obtain R(θ _1b):

(10)

This approach is comparable in some ways to Equation (3), but it has the important advantage of being able to separate out source variability from bed variability. This could be potentially important for exploring, for example, so-called sticky spots beneath glaciers and ice streams. Subglacial friction variability often has surface expression in topography and accumulation variability. Thus in regions of significant topography, this technique could reduce the error introduced by source variability. The disadvantage of using Equation (10) is that like Equation (3) it is sensitive to potentially large uncertainties in the attenuation. Hence in section 4.2 a different approach is explored.

4.2. Estimation of the relative reflection coefficient, rAVA

Here we are interested in measuring the angular dependence of the reflection coefficient, or rAVA (relative AVA). By relinquishing the requirement for an absolute value of the reflection coefficient, uncertainties associated with the unknown attenuation can be substantially reduced. The kernel of the idea is to minimize the uncertainties associated with the attenuation by minimizing the path length differences in the ratio. Starting with Equation (5) and rearranging terms to obtain the reflection coefficient ratio which depends on factors that are measurable or predictable or have relatively small uncertainties,

(11)

Thus, at various offsets, r, the ratio of reflection coefficients at different angles can be measured, with only a small contribution from the unknown attenuation (2d ₁ − d ₂ is small). In practice, one would use small offsets to obtain A ₀ via Equation (6), and the receivers at further offsets would be used to estimate the reflection coefficient ratio, via Equation (11). While individual ratios of reflection coefficients at various angles can be computed using Equation (11), an explicit relationship providing the angular dependence of R is lacking.

To address this, the angular dependence can be obtained by recursively relating each angular ratio to a reference angle. Which reference angle is chosen is arbitrary, but a convenient angle would be the highest angle in the dataset (not counting the datum or data used to obtain A ₀). This recursive referencing can be formally written, starting from the first receiver, j = 1, as

(12)

where the superscript indicates the receiver number, and n ⁻ indicates the previous value of n. In other words, the reflection coefficient at a receiver N, R(θ _1b)^(N), can be determined relative to the reflection coefficient at the arbitrary reference angle, here selected as R(θ _2b ⁽¹⁾), so long as the following holds: the angle of the multiple reflection at N is the same as the angle of the primary reflection at another receiver (shown in Equation (12) as n ⁻). Note that n does not necessarily take on consecutive integer values; in other words, the number of products in Equation (12) is not always N–j.

The method might be best understood by examining a specific example. Consider an experiment with H = 3000 m and a smooth basalt bed with compressional, shear velocities and density respectively of V _p = 5.7 km s⁻¹, V _s = 3.3 km s⁻¹ and ρ = 2.7 g cm⁻³ below an ice layer with the velocity profile (provided by L. Peters from a site near the South Pole) shown in Figure 3. There are 101 receivers spaced at 100 m intervals, yielding angles at the bed ranging from 30⁸ to 90⁸. The zeroth receiver (at r = 0 offset) is used to obtain A ₀. The bed angles of the next five receivers are listed in Table 1.

Fig. 3. Compressional (black) and shear (gray) velocities in the upper (firn) layer used in the example. Velocities are assumed constant below 200 m.

Table 1. Angles at the bed for the first and second bounce. Receivers are spaced at 100 m intervals; source depth 60 m

For the regular receiver spacing here, it can be seen that the bed angle at θ _1b ⁽¹⁾ is the same as the bed angle at θ _2b ⁽²⁾, where the superscript refers to the receiver number. Thus, we can obtain an estimate of R(θ ₁ ^b(2) = 88.1°) relative to R(θ _2b ⁽¹⁾ = 89.5°) by

(13)

The next angle to be estimated is for n = 4 (receiver 4) because the angle of the multiple at receiver 4θ _2b ⁽⁴⁾ is approximately equal to the angle at another receiver (receiver 2) θ _1b ⁽²⁾.Thus, the estimate of R(θ _1b ⁽⁴⁾ = 86.1°) relative to R(θ _2b ⁽¹⁾ = 89.5°) is

(14)

This process can be continued, and, in this case of equally spaced phones, the receivers at which the angles for the multiple match those of the primary for another receiver are multiples of 2, i.e. N = 1, 2, 4, 8, 16… If the value of R(θ _2b ⁽¹⁾) was known, then of course the absolute reflection coefficient would be obtained. However, if R(θ _2b ⁽¹⁾) is not known, one can choose an arbitrary value for it and still obtain a correct estimate of the angular dependence.

In practice, the receiver spacing may not lead to any receiver pairs with the condition θ _1b ⁽²⁾ = θ _2b ⁽ⁿ⁾ precisely met. Nevertheless as long as the angles are reasonably close (the assumption is that R is slowly varying between the two angles), the errors will be modest. If the reflection coefficient is not slowly varying (which might be conservatively assumed if the closest angles differ by many degrees) the reflection coefficient ratios can be interpolated.

For our example of 100 receiver locations (again, the zeroth one is used to estimate A ₀), and starting at receiver No. 1, this procedure will lead to eight estimates (N = 1, 2, 4, 8, 16, 32, 64, 100) of the relative reflection coefficient. However, one can also start at other values, j = 3, 5, 7…, to obtain more estimates of R. In this general case, the method is most succinctly written by starting from the ratio of the reflection coefficient at the angle of interest and moving up in grazing angle (decreasing n)

(15)

While this is superficially the same as Equation (12), it allows for arbitrary values of N, whereas Equation (12) does not. Furthermore, in practice, this approach numerically forces the largest approximation or interpolation near normal incidence where the reflection coefficient should be most slowly varying, hence minimizing error.

To complete the example above, the relative reflection coefficient has been calculated by this method and compared with the exact solution. In Figure 4, the gray crosses represent the numerical solution obtained via Equation (15) with an arbitrary reference of |R(89.5°)| = 0.4. Note that the angular dependence of the rAVA solution closely follows the exact solution with a multiplicative offset. In order to more clearly determine the accuracy of the method, the solution if the normal incidence reflection coefficient was known is also plotted. The excellent agreement indicates that the underlying idea is sound. The small errors in the rAVA estimate (<2% and invisible on the scale of the plot) are associated with the assumption that the R is varying slowly between and when the two angles are not identical. Interpolation could reduce these errors, but seems unnecessary for the synthetic data here.

Fig. 4. Comparison of the exact reflection coefficient of ice over basalt (black line) with rAVA method (gray crosses) without a priori knowledge; an arbitrary reference of R(89.5°) = 0.4 was used. The properly normalized angular dependence (gray circles) agrees very closely with the exact solution, indicating that the recursive referencing procedure is correct.

Though the method does not give the absolute reflection coefficient, the angular dependence of R nevertheless provides a great deal of information about the bed characteristics. In some cases, there are enough clues in the angular dependence to obtain a good estimate of absolute R. For example, in Figure 4, even lacking knowledge of absolute R, it is clear that there is a critical angle associated with compressional waves at 47_°, which (given an estimate of properties at the base of the ice) would yield a good estimate of the compressional speed of the bed (in this case 5700 m s⁻¹). Furthermore, using that value and an appropriate density would give a reasonable value for the normal incidence reflection coefficient and hence yield an estimate of absolute R. In cases where the angular range is limited, such clues may or may not exist.

It should be noted that this approach can also be used to obtain the angular dependence of R _ss, R _ps and R _sp. For the latter two quantities, the normal incidence reflection coefficient is identically zero (hence unobservable), so some other reference angle should be used. It is also noted that while the method gives the magnitude of R, the sign of R can be obtained from the polarity of the bed reflection.

4.3. AVA and rAVA uncertainties

Here the uncertainties of the methods discussed above are briefly explored. In particular, we focus on the unknown attenuation coefficient which in some environments tends to drive the overall uncertainty.

As a numerical example, consider the uncertainties associated with the bounds given by Smith (1997) for ice attenuation. He used α = 0.21 × 10⁻³ m⁻¹, with lower and upper bounds of 0.067 × 10⁻³ and 0.46 × 10⁻³ respectively. Using the environment given above, we can quantitatively examine the contribution of attenuation uncertainties associated with the rAVA approach compared with that of AVA (using Equation (6) in conjunction with Equation (1); see section 4.1). This comparison is shown in Figure 5a, where the uncertainties associated with rAVA (black) are more than an order of magnitude smaller than that associated with AVA (gray). The quantity plotted is the amplitude attenuation factor exp(αD), where D is the effective distance. The greatly reduced uncertainty associated with rAVA is due to the decrease in the D: D = (2d ₁ − d ₂)/2 for rAVA but D = d ₁ for AVA. These distances vs angle are shown in Figure 5b. One can also compare the uncertainties of Equation (4) used by Smith (1997, 2007) and others with rAVA at normal incidence in Figure 5, since Equation (4) has the same uncertainties as AVA at normal incidence.

Fig. 5. (a) Comparison of uncertainties associated with the amplitude attenuation factor exp(αD), where D is the effective distance associated with each method: rAVA (black line) and AVA (gray line). (b) Effective distance associated with rAVA (black line) and AVA (gray line).

It should be noted that the uncertainties of Figure 5a were calculated under the condition that the attenuation is constant with depth. If attenuation is not constant with depth (which is generally the case), the uncertainties will evolve differently with angle than shown here and the uncertainties will be larger. Nevertheless, Figure 5 gives a reasonable first-order picture of the uncertainties for a simple environment.

Inspection of Figure 5 indicates that by using rAVA, the angular dependence of the reflection coefficient is obtained with smaller uncertainties than the absolute reflection coefficient. This should be useful for probing the bed properties inasmuch as the angular dependence may provide many useful clues as to the bed properties. While several strategies could be invoked depending on the nature and quality of the measurements, another possibility would be to use Equation (4) to estimate R(π/2) and then use the rAVA technique to obtain an absolute estimate of R(θ).

It is not the intent here to do a full uncertainty analysis, but note that Equation (15) is the product of ratios derived from the data; this means that errors and uncertainties will propagate. Bias errors will have the greatest impact at low grazing angles (where multiple biased terms are required). Random errors, on the other hand, will be smallest at low angles, since the product of many terms will have an averaging effect. Since rAVA inherently uses ratios, some common causes of bias errors would disappear, including errors in receiver calibration and source amplitude. In order to explore the impact of random errors, a zero-mean Gaussian ‘noise’ with a standard deviation of 0.05 has been added to the reflection ratios (right-hand side of Equation (15)), with the result shown in Figure 6. Besides the averaging effect mentioned above, one would expect larger errors in angular regions where the reflection coefficient ratio is rapidly changing (recalling discussion above on interpolation errors). Simulation shows that for this example, the standard deviation is ∼0.05 at low angles, 0.1 near the peak (critical angle) at 47° and 0.06 near normal incidence. The salient point of the plot is that at this level of noise or uncertainty, the angular dependence of the reflection coefficient still contains a great deal of information. If the random errors are large, their impact can be reduced by averaging over angle or fitting the data to a functional form and interpolating based on the fit.

Fig. 6. Effect of random errors in reflection ratios on rAVA (gray dots). The errors are zero-mean Gaussian distributed with a 0.05 standard deviation. For clarity in presentation, it is assumed here that the normal incident reflection coefficient is known.

As a final comment, it should be noted that the strategy for developing the rAVA technique exploits a ratio A ₁ ²/A ₂ which minimizes the total path length, thus minimizing uncertainties due to the attenuation. At the same time as minimizing attenuation uncertainties, this ratio also tends to reduce uncertainties associated with the spreading factor contained in γ.

5. Conclusions

Analysis of seismic reflection data to obtain estimates of subglacial bed properties requires estimation of source amplitude and attenuation in the ice (among other factors). Generally the attenuation is poorly known and this leads to large uncertainties in the reflection coefficient, i.e. large uncertainties in the bed properties. For example, the uncertainties due to attenuation can be so large that it is not possible to distinguish unconsolidated sediments from crystalline rock.

In this paper, simple methods were derived for extracting the source amplitude and the angular dependence of the bed reflection coefficient minimizing the effect of uncertainties associated with the unknown (or poorly constrained) ice attenuation. The two methods for obtaining source amplitude do not require any knowledge of attenuation in the medium. The method for the angular dependence of the reflection coefficient does require an estimate of the attenuation, but it is much less sensitive to uncertainties in the attenuation than previous methods. These methods are important to the glaciology community because they firstly provide new strategies for processing data to recover conditions at the base of ice sheets that are critical to better numerical modeling. Secondly, the techniques we suggest will aid in the design of new experiments targeted at those environments.