DETERMINATION OF THE EFFECTIVE ELASTIC THICKNESS OF THE CRUST BY SPECTRAL ANALYSIS

by John Moores [990161846]

Submitted: April 24, 2003 Toronto Ontario Instructor: R. Bailey PHY 495S

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

1

TABLE OF CONTENTS

1. 2. 3. 4. 5. 6. 7. 8.

Introduction and Motivation………………………………………. Notation……………………………………………………………………… Theory and Approach…………………………………………………. Implementation………………………………………………………….. Error Considerations………………………………………………….. Results and Discussion………………………………………………. Conclusions………………………………………………………………… References and Bibliography……………………………………..

Page Page Page Page Page Page Page Page

2 4 5 10 14 16 19 20

Appendix A – Graphs and Figures Appendix B – MATLAB Codes

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

2

INTRODUCTION AND MOTIVATION While it is important to understand the stratified division of the Earth into inner and outer core, mantle and lithosphere, one must keep in mind that none of these layers is homogenous. For the surface of the lithosphere this inhomogeneity is as obvious as the differences between plains, mountains and oceans. However, these features are not exclusively the result of internal processes but also interact with the layers below in a feedback loop of sorts. For instance, while tectonic processes deep in the mantle drive the formation of high mountain chains, given enough time the additional bulk will deflect the lithosphere and the entire mountain chain will sink until an hydrostatic equilibrium is reached. This is referred to as isostatic compensation. As such, we are interested not only in the composition and distribution of these layers, but also how they respond to loading. In particular, this paper is concerned with the elastic deformation of the lithosphere in response to the periodic loading of the terrain. Since the layers of the earth near the surface may be approximated as plates, this suggests the definition of an 'effective elastic thickness' as the depth of lithosphere over which the response can be likened to the deformation of a plate under load. However, some caution is warranted. By no means should this be interpreted to be the actual thickness of the lithosphere itself or that there is a discontinuity demarcating a boundary between two layers at this depth (even though we will model the lithosphere in this way). In fact, it is well known from seismic and thermodynamic arguments that the boundary between lithosphere and mantle is somewhat 'blurred' and occurs over several kilometers1. Instead, it is best to conceptualize this effective elastic thickness as the thickness of a perfectly elastic infinite plate composed of lithospheric material located on top of an essentially fluid (at least on geological time scales) mantle. Thus, using the principles of structural/continuum mechanics we can determine this thickness by observing the deformation in the plate due to a known loading pattern. While this may not give us the lithospheric thickness directly, it does allow for the comparison of different regions since it is logical to suppose that if a particular region has a higher effective elastic thickness, it also has a thicker lithosphere. We can learn a great deal about the lithosphere from this thickness variation. For instance, we know from thermodynamics that as the temperature of rock increases its elasticity JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

3

decreases thus we can argue that areas with greater elastic thicknesses are also colder and thus presumably older. It is also reasonable to postulate that mountain building would be easier in a region where the lithosphere is thin and thus more prone to buckling from tectonic stresses. The aligned rows of alternating mountains and valleys of the basin and range region of the Rocky Mountains, is highly suggestive of such a buckled plate – a thin elastic thickness here would tend to reinforce this supposition. Thus, we can gain a number of insights into the nature of the upper mantle and lithosphere by determining this elastic thickness for different regions. In particular, we shall compare several regions in North America including the basin and range region of the Rocky Mountains - thought to be a relatively young region - in addition to the Appalachian Mountains and Canadian Shield - both thought to be very old surface features.

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

4

NOTATION φ

North latitude of a particular site from the Equator (Degrees)

ψ ωx ,ωy ωn ρm , ρc E g λ λchar ωchar t h p q Q

East Longitude of a particular site (Degrees) Frequency of oscillation in eastward and northward directions (Rads/km) Nyquist Frequency (Rads/km) Density of the mantle and crust respectively (kg/m3) Elastic modulus of the lithosphere (N/m2) Gravitational Acceleration (9.81 N/kg) Wavelength (m) Characteristic Wavelength of the lithosphere (m) Characteristic Angular Frequency of the lithosphere (Rads/km) Effective Elastic Thickness of the lithosphere (m) Loading Amplitude of the Topography (m) Horizontal Loading of the lithosphere (N/m2) Vertical Loading of the lithosphere at the surface - i.e. the datum (N/m2) Hydrostatic restoring pressure exerted by the mantle on deflected lithosphere (N/m2) Response Amplitude of the Lithosphere (m) Coherence (dimensionless) Flexural Rigidity of the lithosphere (N·m) Degree of Compensation (dimensionless) Poisson’s Ratio (dimensionless) Cartesian coordinates in spatial domain (km) The number of samples in a spatial or frequency domain The spatial or temporal separation between successive samples The period of the domain, given by N ∆t Matrix indices (dimensionless) Fast Fourier Transform function Fast Fourier Transform function in 2 dimensions Topography matrix in space domain Gravity anomaly matrix in space domain Topography matrix in frequency domain Gravity anomaly matrix in frequency domain Power Spectrum of the topographic map Power Spectrum of the topographic map Cross-Power Spectrum of the topographic and gravity maps The error in the Power Spectrum of the topographic map The error in the Power Spectrum of the topographic map The error in the Cross-Power Spectrum The error in the Coherence

z Coh D C ν x,y N ∆t τ i,j,k fft( ) fft2( ) topo gamap T G PT PG STG δPT δPG δSTG δCoh

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

5

THEORY AND APPROACH The most obvious question is how does one determine this elastic thickness? While we can estimate the mass of mountains, it is not obvious as to how to measure the deflection of the lithosphere. The key however lies in the concept of isostacy. As mentioned in the introduction, this is a form of hydrostatic equilibrium - that is large surface features will sink over time until the weight of adjacent columns is equivalent [Figure 1a]. As a result, large surface features have 'roots.' This mechanism, however, is fundamentally different from elastic support [Figure 1b]. In this case, we do not have an equivalence between adjacent columns of material, in fact, there is an overabundance of mass in the vicinity of features which are elastically supported since the deflection of the lithosphere is less than it would be given Figure 1: Isostatically compensated terrain (a) Elastically supported Terrain (b)

isostatic compensation. As such, it is possible to determine the elastic thickness of the

lithosphere by determining the characteristic wavelength at which the topography shifts from one type of behavior to the other. We may relate this crossover wavelength by considering the differential equation for a plate (Please note that the derivation that follows parallels that found in Turcotte and Schubert2, chapter 3):

D

d 4z d 2z + p = q ( x) − Q ( z ) dx 4 dx 2

[1]

Where p is the horizontal stress in the plate, q is the loading at the theoretical surface of the lithosphere (i.e. at the surface of a hypothetical equipotential surface), Q is the hydrostatic restoring force of the liquid mantle and D is the flexural rigidity given by:

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

6

D=

Et 3 12(1 − ν 2 )

[2]

If we assume that the deflection of the lithospheric plate is ‘filled in’ [ Figure 2 ] by crustal rocks, the buoyancy force exerted by the mantle may be expressed as the weight of mantle displaced subtracting the crustal fill:

Q( z ) = ( ρ m − ρ c ) gz

[3] If we further assume that the horizontal stress in the plate is negligible (p≈0 – may not necessarily be true of a tectonically active region) we derive:

D

d 4z + ( ρ m − ρ c ) gz = q( x) dx 4

[4]

This is the fundamental equation for an elastic lithospheric plate under a continent. We now model the loading pattern q. We choose a sinusoid to accomplish this - recall that Fourier theory states that any arbitrary signal Figure 2: We assume an equipotential surface, thus, at the surface, the deflection appers ‘filled in.’ NOTE: for the purposes of this paper, the lower lithosphere has the same properties as the mantle. [Reprinted from Reference 2]

(such as topography on a continent) can be built up from a sum of sinusoidal frequencies10. Thus let us express the topography height (compared to the datum) for a

particular wavelength as:

h = h0 sin(2π

x

λ

)

[5]

if we assume constant density within the terrain, we may find the loading pattern due to h:

q( x) = ρ c gh0 sin(2π

x

λ

)

JOHN MOORES

[6]

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

7

Furthermore, since the loading is sinusoidal we know from experience that the response of the elastic plate will be also, thus:

z ( x) = z 0 sin(2π

x

λ

)

[7]

Substituting [7] into [4] we can derive an expression for the response of the lithosphere to a given terrain load:

z0 =

h0

ρn D ⎛ 2π ⎞ −1+ ⎜ ⎟ ρc ρc g ⎝ λ ⎠

[8]

4

Note that as λ→∞ we have:

z 0 (λ = ∞ ) =

ρ c h0 ρn − ρc

[9]

If this equation is examined closely it can be shown that cross multiplying the fraction and multiplying both sides by g yields q(x) = Q(x). That is to say the additional weight of the terrain is equal to the weight of mantle displaced. Recall that this is the definition of an isostatically compensated surface feature. As such we may define the degree of compensation as:

h0 4

ρm D ⎛ 2π ⎞ −1+ ⎜ ⎟ z 0 (λ ) ρc ρc g ⎝ λ ⎠ ρm − ρc = C (λ ) = = 4 ρ c h0 z 0 (λ = ∞ ) D ⎛ 2π ⎞ ρm − ρc + ⎜ ⎟ ρm − ρc g⎝ λ ⎠

[ 10 ]

It can be seen from graphing members from this family of equations [Figure 3] that there is a rapid cross-over from C ≈ 0 to C ≈ 1. Since this represents a shift from elastically supported to isostatically compensated behavior, we define the characteristic wavelength as being the wavelength corresponding to C(λchar) = 0.5. Substituting these into equation 10 we derive:

⎛ 2π ⎜⎜ ⎝ λchar

4

⎞ g (ρ m − ρc ) ⎟⎟ = D ⎠

JOHN MOORES

[ 11 ]

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

8

Finally, substituting for the flexural stiffness using equation 2 and rearranging gives us an expression for the thickness of the lithosphere in terms of this characteristic wavelength:

t3 =

12 g (1 − ν 2 )( ρ m − ρ c ) ⎛ λchar ⎞ ⎟ ⎜ E ⎝ 2π ⎠

4

[ 12 ]

For this calculation we will use some typical values for the lithosphere which we shall assume is both homogenous and isotropic. These values are given in table 1: TABLE 1: Typical Values for the Lithosphere and Mantle Elastic Modulus2 70x109 N/m2 Mantle Density7 3300 kg/m3 7 Crust Density 2800 kg/m3 7 Lithosphere Poission ratio 0.25

Thus, the original problem of directly finding the thickness of the lithosphere has been reduced to determining the characteristic wavelength that divides compensated from uncompensated behavior. To find this value, it is important to remember that if the terrain is elastically supported there is a surplus of mass in columns which include the raised topography compared to the surrounding countryside. Conversely, for isostatically compensated terrain, all vertical columns have the same mass. Thus we would expect that for small-wavelength topography we would observe a positive gravity anomaly with the same wavelength but no anomaly when considering compensated long-wavelength topography. We must be careful, however, since there is an implicit assumption in this view of the theory. We have assumed that we are measuring the gravity anomaly at a specific constant altitude. This is commonly referred to as the ‘free-air’ gravity anomaly and is the ideal dataset to use. Unfortunately, it is much more difficult to obtain free-air gravity readings – which require an extremely stable platform at altitude1 - then it is to simply walk over the terrain and take stationary readings every few miles. This latter scheme of measurement obtains what is called the bouguer gravity anomaly and also includes an implicit assumption. Specifically, the bouguer gravity anomaly data is not raw data, but has been processed to compensate for the changing altitude of the observer.

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

9

Typically, this is performed by subtracting out the estimated contribution from the extra mass between the observer and the datum (i.e. the mountain they are standing on). Thus, the bouguer gravity anomaly assumes that all terrain is elastically supported. As such, at short wavelengths of topography we observe no corresponding gravity anomaly. However, at large wavelengths, the terrain is compensated and thus the correction to the raw data produces an artificial gravity low at the same wavelengths as the compensated terrain. This is sometimes referred to as the ‘hollow mountain’ problem – since the gravity anomaly produced by large, compensated mountain chains is lower then would be expected if they were elastically supported. This paper takes advantage of this artifact of the data. There is some debate between supporters of either bouguer or free-air anomalies as to which method is the most accurate. Both have problems. The free-air anomaly relies heavily on Figure 3: Degree of Compensation family of curves. The degree of Compensation is plotted against wavelength and frequency (1/λ). NOTE: crossover occurs at C = 0.5.

accurate altitude keeping and positional awareness while also requiring extremely stable inertial platforms for the gravimeters on

notoriously unstable (when compared to the requirements) aircraft. The bouguer data, while more readily available then free-air data and potentially more accurate, require data for large sections of terrain to detect the long-wavelength artifact. Thus, it may be impossible to obtain accurate local thickness estimates. This is only intended to give a flavor to the situation – the resolution of this debate is beyond the scope of this paper. Instead, we will use the bouguer data exclusively. JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

10

In either case, we proceed by performing a spectral analysis on the terrain and gravity maps in what is known as the Coherence Method5,6. Once these have been transformed into the frequency domain we compute the power spectrum of both maps as well as the cross power. These are defined as follows:

PT ≡ (1/n) < Tij* · Tij >

[13-15]

PG ≡ (1/n) < Gij* · Gij > STG ≡ (1/n) < Tij* · Gij > Note that the individual power spectrums of the maps will be real (since they are the product of individual elements of transformed maps multiplied by their complex conjugates) while the cross-power, in general, is complex. From this we may define the Coherence as:

Coh =

S TG

2

[ 16 ]

PT ·PG

It can be shown that the coherence is 1 when the frequencies in T and G are correlated and 0 when they are uncorrelated. Thus we expect that for bouguer gravity the Coherence will be perfect at low frequencies (long-wavelength) and fall off to zero at the characteristic frequency. For the free-air anomaly, we would expect that the Coherence will be zero at low frequencies and jump to perfect coherence at the characteristic frequency. For perfect datasets, we would expect both estimations of the characteristic frequency to be identical, however, in practice this does not occur4.

IMPLEMENTATION Using the MATLAB programming environment, code was written to calculate the value of the coherence over the averaged power spectrums. Three geographic regions were investigated – the Basin and Range, the Appalachians and the Canadian Shield. The source data for the topography in the form of USGS tiles is displayed in Appendix A.1 with the approximate extracted regions highlighted with colored overlays. The actual extracted regions are displayed in Appendix A.2. This source data was obtained at http://edcdaac.usgs.gov/gtopo30/gtopo30.html and was degraded from its initial 30 arc second grid spacing to 2.5 arc minutes8 in order to match with the bouguer gravity map.

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

11

Likewise, the bouguer gravity anomaly map is displayed in Appendix A.3 with the specific extracted regions displayed in Appendix A.4. This was obtained from the binary raster image available from NOAA at http://www.ngdc.noaa.gov/seg/fliers/se-2004.shtml . Two notes of caution are required in the context of this dataset. First, in order to better display the map, the image was created by reversing the bit ordering (given format is 16 bit LSB) to MSB/big endian in order to exaggerate the differences between elements. For the purposes of coding, however, the original little endian bit ordering was employed. Secondly, it was found (though a significant effort) that there is an error in the header file9. The data is claimed to begin at 0.75°W and end at 180°W but instead was found to be deviated to the east by close to 10°, thus the original gravity maps of basin and range were actually being produced from data recorded in the middle of the pacific. It is unknown if this error is unique to the raster file or if it is replicated in the other data formats in the NOAAprovided package. It is important to discuss the criteria for the selection of patch dimensions. Table 2 enumerates the locations and dimensions of the selected patches:

feature Latitude limits Longitude limits Maximum λchar

TABLE 2 Patch Sizes and Locations Patch Name Basin and Range Canadian Shield 37°N – 42°N 45°N – 65°N 113°W – 118°W 94°W – 99°W 557km 2230km

Appalachian 45°N -60°N 61°W – 76°W 1670km

The basic criteria are these: the patches must be as large as possible since we are most interested in correlating long-wavelength topography. However, at the same time we should attempt to limit the number of different geological regions included in a particular patch - that is, we assume that the characteristic wavelength is a property of a particular geological region and is different for adjacent regions. Thus we have an upper limit on the size: for instance the basin and range patch cannot be larger then the basin and range geological region. This led to some peculiarly sized patches. In particular, the Canadian shield patch depicted in Appendix A.1 which attempts to reduce the minimum observable frequency (maximize the observable wavelength) at the expense of the number of elements to be averaged since the slice is very thin in longitude. Furthermore, since there is no elevation data for water, patches containing ocean should be avoided if at all possible. In the case of the Appalachian patch we have ‘cheated’ a little, and it is expected that this will introduce some error, however, it was JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

12

thought that the selected patch would give a better picture of the region then any individual smaller slice and that since the surface area of ocean included is small this effect could be minimized. A useful MATLAB toolbox, gtopo30, was available to read in and process the data tiles for the USGS data set, however, since the NOAA data was significantly older it was necessary to write a separate program to read in and order the data. This file – ggrav.m - is provided along with the rest of the code in Appendix B. Once this data was successfully introduced into the MATLAB environment, a bit more preprocessing was required. In order to eliminate DC power problems and broadening at low frequencies, the mean was removed and the resulting maps were windowed. Three types of windows were attempted – Welch, Hanning, and Bartlett11. These are displayed for comparison in Appendix A.6. While all three were found to be effective, the Welch window appeared to perform marginally better then its counterparts and thus was used on data from all three regions. Following this, the maps were transformed into the frequency domain using the built-in Matlab 2-dimensional fast fourier transform function (fft2) before being reduced. A bit of explanation is required regarding the way in which transformed matrices and vectors are stored in the frequency domain. The transform of a discretely sampled signal has the same number of elements as the original signal. For a vector input, we may obtain the frequency of any frequency domain sample according to:

ωk =

2π 2π k= k τ N∆t

[ 17 ]

Where k is the index of the vector and τ is the spatial or temporal period of the original time or spatial signal – that is the number of samples (N) multiplied by the temporal or spatial spacing between each sample (∆t). From discrete fourier theory, we know that the maximum frequency which we can extract from a time signal has a wavelength of two divisions, that is:

ωn =

2π 2∆t

[ 18 ]

This maximum frequency is called the nyquist frequency. Note that we obtain the nyquist frequency when k = N/2 – this begs the question: what is stored in the rest of the vector? As it turns out, this stores the spectrum from the negative nyquist back to zero1 [ Figure 4 ].

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

13

This bookkeeping (since negative frequencies are no less physical) means that we have stored frequencies in the range of – ωn to + ωn even though the negative portion of the spectrum is redundant. In fact, for visualization purposes, the function fftshift in MATLAB will reorganize the vector such that the transform is symmetric about the zero frequency instead of the nyquist. This situation is more complex in two dimensions [ Figure 5 ], but remains analogous. Here we have four quadrants in which the first and third are copies, as are the second and fourth, however, each set contains independent information. Thus it is necessary to perform a superposition to ensure that all the pertinent data is extracted. For the purposes of reducing this data, we superimpose the first and second quadrants (the elements of the second are reversed in Figure 4: Array ordering in a 1-Dimensional Fourier Transform

ωx, but not in ωy). This ensures that the averaging routine which follows

(concentric circles) has all the independent information. First, we obtain the power spectra by employing equations 13-15.

Figure 5: Array ordering of a 2-Dimensional Fourier Transform showing concentric averaging circles and the effect of fftshift. JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

14

Next we average values along lines of constant frequency. Since we have assumed that the lithosphere is isotropic this was accomplished by averaging in concentric circles. In order to keep the routine as general as possible, an elemental averaging was performed with ωx scaled with respect to the ωy by the spatial period at the middle latitude in the range. Since both maps provided each have a geographic projection this means that while the absolute size of the divisions in the latitude, φ, do not change, the size of the divisions in the longitude get smaller as we approach the poles. As such, scaling by the average is at best an approximation. This leaves only the plotting of the coherence values. For all three zones it was decided to plot averages over 1,2 and 3-elements. This refers to the number of elements in y between concentric circles in the averaging routine. The resulting coherence patterns have been plotted in Appendix A.4 on a semi-log scale vs. the corresponding angular frequency.

ERROR CONSIDERATIONS Before the results of the analysis are discussed, a digression regarding error considerations is required. First, qualitative concerns will be discussed which will be followed by an attempt to quantify the statistical error involved in the coherence calculations assuming that the data up to that point have no intrinsic measuring error – that is they are assumed to be perfect samples of a stochastic process. The first consideration has to be for the source data. Since we are dealing with distributions, we have two potential errors – measurement errors and positioning errors. Generally speaking, the topographic map appears to be more reliable since it is taken from relatively few sources, has a greater resolution and is significantly more recent then the gravity anomaly map8,9. Furthermore, measuring elevation, especially on a large scale, is a far simpler task then measuring the gravity anomaly. Additionally, we know that there is at least one systematic error in the gravity raster file which was detected – the entire map was shifted by close to 10° to the East. The correct axes were only obtained by referencing easily identifiable locations on the gravity map (even these were only available since Mexico was a ‘no data’ region and the map could be calibrated by

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

15

locating Yuma and Cancun) Despite these obvious flaws, for the purposes of this paper, we shall consider both these errors (positional and measurement) to be zero for both topographic and gravity anomaly maps. Lastly, this is not to say that the maps are not noisy. In fact, the bouguer gravity anomaly mapping technique was originally devised to search out density variations in the crust for the purposes of oil, gas and mineral exploration1. Since we are interested in long wavelength features, these density variations are simply high-frequency noise and can be ignored. The second consideration results from the projection used and from the spherical nature of the earth’s surface. For the purposes of this project, it is assumed that the selected regions are perfectly flat . This is known to be false - especially for the larger plates which can measure in excess of 2000km on a side - but may not be a bad approximation. However, the gridding of the data is of much greater concern. The longitudinal period was estimated by using the value at the median latitude and assuming that the regions were rectangular. For certain plates this is highly questionable, especially in the case of the Canadian Shield Patch. Here the period in y is 20° (45°N to 65°N) and in x it is 5°. The scaling of the longitude domain is proportional to the cosine of the north latitude, thus the employed scale factor (at 55°N) is 0.5736, however the actual variation from the bottom to the top of the patch varies from the high of 0.7071 (at 45°N) to the low of 0.4226 (at 65°) – a change of just over 40%. Naturally this is the worst case – these figures are much lower for the basin and range patch since it is much smaller and located closer to the equator. This needs to be taken into account when comparing the reliability of the result obtained from different patches. In the special case of the Appalachian patch we must also be mindful of the inclusion of ocean surface -as was alluded to in the last section. Even if we neglect all these errors, we must still consider the statistical error in getting from the various power spectra to the coherence. From basic statistics we may define these errors (characterized by calculating the standard deviation) as:

δPT = δPG =

1 2NT 1 2N G

PT

[ 19 ]

PG

[ 20 ]

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

16

δS TG =

1 2 N TG

PT PG

[ 21 ]

Thus by the SSE method we have, for the standard deviation in the Coherence:

δCoh =

∂Coh ∂Coh ∂Coh δS TG + δPT + δPG ∂S TG ∂PT ∂PG =

2 S TG

1

PT PG

2 N TG

PT PG +

− S TG 2

PT PG

2

1 2NT

PT +

− S TG PT PG

2 2

1 2N G

[ 22 ]

PG

If desired, the code could be extended to take into account this statistical error, however, it is provided here purely for illustrative value.

RESULTS AND DISCUSSION The coherence distribution is given in Appendix A.5. From these graphs, the characteristic frequencies are obtained for the fall-off where the Coherence is 0.5. This is summarized in table 3 in which we have assumed the values from Table 1 in order to calculate the elastic thickness: TABLE 3: Summary of Characteristic Attributes of the three Patch Name Attribute Basin and Range Canadian Shield ωchar 5-8.5x10-3 rads/km 6.5x10-2 rads/km λchar 96.7km 1300km – 740km Elastic thickness 3.5 km 110km – 53km

Patches Appalachian 0.9-1.5x10-2 rads/km 700km – 420km 49km – 25km

The first thing to notice about this data is the roll-off properties exhibited by the figures in Appendix A.5. The curve for the Basin and Range looks a great deal like the expected theoretical form of Figure 3. We have a clear plateau and a sharp fall-off. In this case the most likely curve is the 2-point average, however, all three curves are very similar at the cross-over coherence.

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

17

The plateau feature is much less pronounced on the graphs for the Appalachians and the Canadian Shield. This is also expected since these regions were anticipated to produce long characteristic wavelengths due to large thicknesses. As a comparison, note that characteristic wavelength for Basin and Range of 96.7km (table three) is only 17% of the maximum value (table 2), given the patch size while this figure is 57% for the Canadian Shield. Unfortunately, since the geographic domains of the particular regions are limited it is impossible to increase the patch size to encompass more of the roll-off curve. One particularly puzzling aspect of the curves is the presence of a spike in the data just after the fall-off. This is particularly pronounced in the basin and range data and does not seem to be eliminated by increasing the averaging width. This suggests that this might be physical. If this is the case, it is possible that we are seeing local variations within the region. For instance, if, locally, within the basin and range zone there was a region with a substantially thinner section of elastic lithosphere (for instance due to volcanism or a hotspot) we could have compensated geography in a section of the spatial domain at very small wavelengths – this would give rise to an improved coherence at these higher frequencies. This could also be accomplished without requiring a change in lithosphere elastic thickness if the lithosphere were to locally suffer a failure (i.e. to break from loading) and thus lose much of its flexural rigidity. This is not an altogether impossible situation – note that the thickness of the lithosphere in the basin and range region is comparable to the height of the topography suggesting that there may be buckling failures at specific locations within the region which allow higher-frequency topography to be compensated. Despite this speculation, it is equally possible that this is simply an anomaly of the data sets used and has no physical significance. It is interesting to note that the ordering in terms of thickness is what was expected with the Basin and Range patch exhibiting the smallest elastic thickness and the Canadian Shield exhibiting the largest. These thicknesses are also compatible with results obtained by Bechtel et al3 using a coherence method on bouguer gravity data. That group obtained an average value of 5km for the basin and range and a minimum of 4km which agrees well with the value obtained here. Furthermore, for the Canadian Shield they found a maximum thickness of 128km and an average thickness of 82km while for the Appalachian region a thickness of 39km was obtained. These also compare well with the values in Table 3. While this is encouraging, it is important to recognize that Bechtel et al possessed bathymetric data for Hudson’s bay and included this in their average over the Canadian Shield. Furthermore, since the code provided here is intended to deal with rectangular JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

18

patches parallel to the latitude and longitude axes the patch used for the ‘Appalachians’ actually extends over three geological regions, namely the Canadian Shield, Appalachians, and Grenville regions [ figure 6 ]. However, Bechtel’s Appalachian average is composed of only the latter two regions and also uses Bathymetric data for the Gulf of St. Lawrence to obtain a patch of sufficient size. The largest discrepancies appear to occur for the Canadian Shield and Appalacian region. This is not surprising – as discussed in the section on error we have more reason to distrust our values from these patches then the Basin and Range region. As well, it should also be noted that Bechtel’s data has been re-projected on a Lambert conformal projection and avoids the error of Figure 6: Geological Map of North America [Reprinted from Reference 3 ]

skewed patches at northern latitudes to which the code described in this paper is

subject. As such, the method outlined here could be made more precise by re-gridding the data to a Cartesian plane – this is the simplest extension of the code. It would also be beneficial to examine free-air gravity anomaly maps. Should these be available the current code would be able to analyze the coherence with the simple addition of a utility to read in the data for processing. Not only would this give us a lower bound for λchar but might be obtained for smaller patches. This is due to the roll-off properties of the free-air data which occurs at lower wavelengths then for the bouguer data4 (of course, given perfect data these would be identical, however, experiment suggests that the bouguer gives an upper bound while the free-air method gives a lower bound) and the fact that we have more highfrequency, small-wavelength data then low-frequency data suggests that the roll-off might be more pronounced then the bouguer graphs in regions with greater elastic thicknesses.

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

19

CONCLUSIONS Using a coherence method, values were obtained for the elastic thickness of the lithosphere in the Basin and Range, Canadian Shield and Appalachian regions of North America under the assumption of an isotropic lithosphere. The values obtained for the characteristic wavelength were 3.5km, 110-53km and 49-25km, respectively, in each region. These values agree well with the data obtained by Bechtel et al especially considering the concerns about the method which we have used at higher latitudes on larger patches. Furthermore, they support the hypothesis that the lithospheric thickness is related to the age of the surface features observed. This code could easily be extended and the accuracy of the results improved substantially by re-gridding the data on a Cartesian projection and by comparing the characteristic wavelength to that obtained by using a free-air method.

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

20

REFERENCES [ 1 ] Bailey, Richard. (2003) PHY 495 Course Notes. University of Toronto. [ 2 ] Turcotte,D. L. and G. Schubert (2002) Geodynamics. Second Edition. Cambridge University Press. [ 3 ] Bechtel T. D. et al (1990) “Variations in effective elastic thickness of the North American Lithosphere.” Published in Nature vol 343 pp 636-638. [ 4 ] McKenzie, D. and Fairhead D (1997). “Estimates of the effective elastic thickness of the continental Lithosphere from bouguer and free-air gravity anomalies.” JGR-Solid Earth paper 97JB02481, vol 102 n°12, pp.27523. [ 5 ] Forsyth, D.W. (1985) “Subsurface Loading and estimates of the flexural rigidity of continental lithosphere.” Published in the Journal of Geophysical Research, vol 90, n°B14 pp.12623-12632. [ 6 ] Macario, A. et al (1995) “On the robustness of elastic thickness estimates obtained using the coherence method.” Published in the Journal of Geophysical Research, vol 100, n°B8 pp.15163-15172. [ 7 ] Sandwell, D.T. (2002) Gravity/Topography Transfer Function and Isostatic Geoid Anomalies. Notes. [ 8 ] United States Geological Service (1998) gtopo30 supporting documentation. [online] http://edcdaac.usgs.gov/gtopo30/gtopo30.html. [ 9 ] National Oceanic and Atmospheric Administration (1989) Bouguer gravity anomaly supporting documentation. [online] http://www.ngdc.noaa.gov/seg/fliers/se2004.shtml . [ 10 ] Haberman, Richard (1998) Elementary Applied Partial Differential Equations with fourier series and boundary value problems. Prentice Hall Publishing [ 11 ] Bourke, Paul (1998) Windows. Swinburne centre for astrophysics and supercomputing. [online] http://astronomy.swin.edu.au/~pbourke/analysis/windows/.

BIBLIOGRAPHY [ 1 ] Serway, R.A. (1996) Physics for Scientists and Engineers. Volume 1. 4th Edition. Saunders College Publishing. [ 2 ] Van Loan, C.F. (2000) Introduction to Scientific Computing: A Matrix-Vector Approach Using Matlab. 2nd Edition. Prentice Hall Publishing. [ 3 ] Nyborg, D. (1998) The Fast Fourier Transform and its use in Spectral Analysis of Digital Audio. Published Online by www.Mathtools.net. [ 4 ] Carrol, B.W. (1996) An Introduction to Modern Astrophysics. Addison Wesley Publishing Inc. [ 5 ] French, A.P. (1971) Vibrations and Waves. Published by WW Norton & Company. [ 6 ] Karner, G.D. and Watts, A.B. (1983) “Gravity Anomalies and Flexure of the Lithosphere at Mountain Ranges.” In the Journal of Geophysical Research, vol 88 n°B12. pp. 10449-10477. [ 7 ] Uncredited Reading. Chapter 6, sections 6.6 to 6.13.

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

21

APPENDIX A

CONTENTS: GRAPHICAL Topographical Maps (USGS Tiles) Topographical Maps (Generated Patches) Main Gravity Map (NOAA Raster Image) Gravity Maps (Generated Tiles) Coherence vs. Angular Frequency Windows

A.1 Topographical Maps (USGS TILES)

Basin

2

& Range Region

Canadian Shield Region

1

Appalachian Region

3 Counterclockwise from Above: 1. USGS Tile w100n90

2.USGS Tile w140n90

3.USGS Tile w140n40

Designation refers to top left hand corner of the tile. NOTE: while the Appalachian and Canadian Shield regions derive solely from w100n90, Basin and Range stretches over tiles w140n90 and w140n40. Approximate areas as shown.

A.2 Topographical Maps (GENERATED PATCHES)

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

24

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

25

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

26

A.3 Main Gravity Map (NOAA RASTER IMAGE)

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

27

A.4 Gravity Maps (Generated Tiles)

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

28

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

29

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

30

A.5 Coherence vs. Angular Frequency

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

31

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

32

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

33

A.6 Windows

JOHN MOORES

DETERMINATION OF THE ELASTIC THICKNESS OF THE CRUST USING SPECTRAL ANALYSIS ( 2003 )

PAGE

34

APPENDIX B

CONTENTS: CODES Matlab Code – Core program Matlab Code – Supporting Functions (general) Matlab Code – Supporting Functions (region-specific)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % PHY 495S ASSIGNMENT #2 % % % % The Elastic Thickness of the Earth’s Crust % % % % by: John Moores (c) February - April 2003 % % Supervisor: R. Bailey % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % MAIN PROGRAM CODE % % CANADIAN SHIELD ZONE % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Lowercase letters indicate space domain quantities % Uppercase letters indicate frequency domain quantities % % The code is identical for all zones except different loaders are used % (APtopo for the Appalachians and BRtopo for Basin and Range) clear all CStopo

%% Script which loads the two matrices into the Matlab Environment %% map = topographical map , gamap = gravity anomaly map

xRange = lonlim(2)-lonlim(1); yRange = latlim(2)-latlim(2); nx = length(map(1,:)); ny = length(map(:,1)); nxq = floor(nx/2)+1; nyq = floor(ny/2)+1;

%% nyquist elements in x and y

%GLAZIER’S SECTION – WINDOW PREPARATION for i = 1:nx %% welch window processing for j = 1:ny welch(j,i) = (1-(((i-1)-(nx-1)/2)/((nx-1)/2))^2 ) *(1-(((j-1)-(ny-1)/2)/((ny-1)/2))^2 ); end end for i = 1:nx %% hanning window for j = 1:ny hanning(j,i) = sin((i-1)*pi/(nx-1))*sin((j-1)*pi/(ny-1)); end end for i = 1:nx for j = 1:ny

%% bartlett window

if(j