Geophysical Exploration Techniques

Materi Kuliah Metode Geofisika Eksplorasi

Fisika-Universitas Fisika-Univers itas Haluoleo 2008

Geophysics 224 - Geophysical exploration techniques Winter 2006 Instructor

Professor Martyn Unsworth Offices : Physics 545B and Physics 605 e-mail : [email protected] Tel : 492-3041

Teaching Assistant

Ersan Turkoglu Office : Physics P605 e-mail : [email protected]

Course Schedule

Classes Lab Mid-term exam Final exam

Tuesday and Thursday Wednesdays th Tuesday February 14 2006 th Thursday, April 27 2006

08:00 to 09:20 14:00 to 17:00 08:00 to 09:20 09:00 to 12:00

P631 P145 P631 P631

Objectives of class

To learn the basic principles of three geophysical exploration techniques. This will include both the physical basis  of the method and data interpretation. Simple mathematics will be used where needed, but the emphasis will be on understanding concepts and methods. Where appropriate, PC based software will be used used in the lectures and the labs. Course outline

A. Introduction B. Gravity exploration C. Seismic exploration D. Magnetic exploration

~ 2 weeks ~ 6 weeks ~ 2 weeks

Office hours

If you have questions about class material or assignments, I will happy to arrange a convenient time to meet with you. To do this, please see me after class, e-mail me or drop by my office. Assignments

There will be a total of four four assignments in this class. These will be due at 5 pm on the appropriate day, and should be handed in at my office (Physics 545B) Recommended text rd   An introduction to geophysical exploration, by P. Keary, M. Brooks and I. Hill, 3 Edition, Blackwell, 2002. ISBN 0-632-04929-4

Course weighting Labs Assignments Mid-term exam Final exam

20% 10% 25% 45%

Statement on academic integrity

"The University of Alberta is committed to the highest standards of academic integrity and honesty. Students are expected to be familiar with these standards regarding academic honesty and to uphold the policies of the University in this respect. Students are particularly urged to familiarize themselves with the provisions of the Code of Student Behaviour (online at www.ualberta.ca/secretariat/appeals.htm) and avoid any  behaviour which could potentially result in suspicions of cheating, plagiarism, misrepresentation of facts and or participation in an offence. Academic dishonesty is a serious offence and can result in suspension or expulsion from the University." (GFC 29 SEP 2003)

Grading in Undergraduate Courses

Description Letter Grade Grade Point Value  _______________________  __________________________________ _____________________ __________ Excellent A+ 4.0 A 4.0 A3.7  _______________________  __________________________________ _____________________ __________ Good B+ 3.3 B 3.0 B2.7  _______________________  __________________________________ _____________________ __________ Satisfactory C+ 2.3 C 2.0 C1.7  _______________________  __________________________________ _____________________ __________ Poor D+ 1.3 Minimal Pass D 1.0  _______________________  __________________________________ _____________________ __________ Failure F 0.0

Policy about course outlines can be found in Section 23.4(2) of the University Calendar (CGF 29 SEP 2003)

MJU 2006

Geophysics 224 - Geophysical exploration techniques A: Introduction What is geophysical exploration? ●imaging underground structures from surface measurements of artificial or natural signals that travel through, or are generated genera ted by, the Earth. Why? ● oil and gas exploration ● mineral exploration ● hydrogeology ● monitoring contamination and remediation ● tectonic studies ● earthquake hazard studies Some basic principles

All geophysical methods remotely sense a material property of the Earth ( e.g. seismic velocity, rock density, electrical resistivity, magnetization etc)

●

Knowledge of these material properties must then be interpreted  to determine what rock type is present. Well log information is very important in this task.

●

Geophysical methods can be divided into active and passive techniques.

●

In an active technique, it is necessary to generate a signal  ( e.g. in seismic reflection surveying, sound waves are generated with an explosion). In a passive technique a naturally occurring signal is detected (e.g.  the pull of gravity of a buried object) Geophysical and geological studies complement one another. Geologists are more effective with a basic knowledge of what geophysics can and cannot resolve. (Similarly, many geophysicists would benefit from a basic knowledge of geology).

●

Geophysical imaging does not always give a unique answer! Additional information is often needed to discriminate between possible solutions ( e.g.  other geophysical surveys, knowledge of local geology, well log information in the study area).

●

Please do not be intimidated by equations! Mathematics will be used in this class where needed, and I hope that this will provide a review of math classes you have taken in the recent past. I do not expect students to memorize equations. My expectation is that students will be able to perform simple rearrangement of equations, and use a calculator to evaluate an equation for a given set of values.

●

I will appreciate feedback to keep the mathematics at an appropriate level.

200 km diameter impact crater in the Yucatan Peninsula of Mexico, imaged with gravity data

Paleo-meanders Paleo-meanders imaged in 3-D seismic data, Gulf of Thailand (AAPG memoir 42)

Plate tectonics in action. Seismic reflection reflection and electromagnetic images of India underthusting the southern edge of the Tibetan Plateau (Unsworth et al., 2005) Unsworth, M.J. A.G. Jones, W. Wei, G Marquis, S. Gokarn, J. Spratt, Crustal rheology of the Himalaya and Southern Tibet inferred from magnetotelluric data, 438, 78-81, doi:10.1038/nature04154, Nature 2005.

B. GRAVITY EXPLORATION B1.1 Newtonian gravitation

Consider two point masses that are a distance r   apart. Newton’s theory of gravitation  predicts that they will attract each other with a force F  force F  that  that is given by:

 F 

=

Gm1 m2 r 2

The quantity G is called the gravitational constant (or “big G”) but is actually a small number. Newton deduced this equation from observing the motion of planets and moons in the solar system. The units are as follows:  F  Newton (N) m kg metres (m) r  2 -2 G = 6.67 x 10 -11 N m kg It can be shown that if one of the masses is finite in size ( e.g. a e.g. a planet or the sun), then  F  due  due to will be the same as if all the mass were concentrated at the centre.

If m1 = m2 = 1,000,000 kg (1000 tonnes) then the variation of  F  as r increases will show the “inverse square law”. law”.

Consider the mass m2. Newton’s Third Law of Motion  predicts that this mass will accelerate with an acceleration a, where:  F = m2 a Rearranging this equation gives

a

 F  =

m2

 Now we know the value of  F  from the first equation so

a  g  =

=

Gm1m2 1 r 2

m2

=

Gm1 r 2

where g is called the gravitational acceleration . Let us now consider that m1 is the Earth and m2 is a small object that we are going to drop. This equation tells us that the acceleration does not depend on the mass of the object being dropped. This was proved by Galileo who allegedly dropped masses masses from the leaning tower of Pisa in Italy. This result says that a small mass and a large mass will fall with the same acceleration.

B1.2 Density of rocks and minerals The Greek letter rho ( ρ) is used to represent density. Rocks and minerals found on Earth have densities that range from 1000-7000 kg m -3. Often densities are quoted in g cm -3. To convert, remember that 1 g cm -3 = 1000 kg m -3. Pure minerals can exhibit a high density since the atoms are closely packed together.

Magnetite Pyrite Galena

ρ = 4.90-5.20 g 4.90-5.20  g cm-3 4.90-5.20  g cm-3 ρ = 4.90-5.20 g ρ = 7.40-7.60 g 7.40-7.60  g cm-3

Sedimentary rocks generally have lower densities since the atoms are not as closely packed together and pore space is filled with lower density materi als.

Water   Alluvium Shale  Limestone  Dolomite

1.00-1.05  g cm-3 ρ = 1.00-1.05 g 1.96-2.00  g cm-3 ρ = 1.96-2.00 g ρ = 2.00-2.70 g 2.00-2.70  g cm-3 2.60-2.80  g cm-3 ρ = 2.60-2.80 g ρ = 2.28-2.90 g 2.28-2.90  g cm-3

The range of density values reflects the degree of weathering and the porosity. Why? The composition is secondary, but explains why dolomite is more dense than limestone. Most sedimentary rocks increase in density with depth, owing to increased compaction. This reduces the pore space available for low density materials such as air, water and hydrocarbons. Igneous rocks are generally more dense owing to minimal porosity.

Granite  Basalt 

2.50-2.70  g cm-3 ρ = 2.50-2.70 g ρ = 2.70-3.20 g 2.70-3.20  g cm-3

The density depends primarily on the rock composition. Mafic rocks are generally

more dense than felsic rocks owing to increased proportion of heavier elements such as Fe and Mg.

B1.3 Units for gravity measurements (milligals)

The gravitational acceleration at the Earths surface is 9.8 ms-2. Subsurface variations in rock density produce very small changes in this value, so it is more convenient to use a smaller unit. 9.8 ms-2

= 980 cm s-2 = 980 gal (after Galileo) = 980,000 milligals

1 milligal

= 10-5 ms-2

B1.4 Approximate calculation to estimate the size of gravity anomalies

Consider a spherical ore body with density  ρ and radius a that is buried at a depth  z below  z  below the surface. The extra pull of gravity will be greatest at a point P directly above the ore body. Remember that the pull of gravity ( g) for a sphere is the same as if all the mass were concentrated at the centre. Total mass of ore body = volume volume  x density = Excess mass , mE

4 3

3

π a  ρ 

= mass of ore body – mass of rock that was already there =

=

4 3 4 3

3

π a  ρ   -

π a

3

( ρ 

−

4 3

3

π a  ρ 0

 ρ 0 )

Above ore body, the change in gravity due to the ore body (gravity anomaly) is given by

Δg

=

Gm E   z 2

=

4Gπ a 3 ( ρ 

−

ρ 0 )

3 z 2

 Now guess some values a  z   ρ  ρ0

= 30 m = 40 m = 4000 kg m -3 = 2000 kg m -3

Careful use of a calculator gives Δg = _____ mgal

 MJU 2006 

Geophysics 224 B2. Gravity anomalies of some simple shapes B2.1 Buried sphere

Gravity measurements are made on a surface profile across a buried sphere. The sphere has an excess mass  M S   x, the vertical S  and the centre is at a depth  z . At a distance  x, component of g is given by

 g  z  =

GM S  z  3

( x

2

2

+ z 

)2

This curve is drawn below for a sphere with: Radius, a = 50 m Depth, z Depth, z -3 Density contrast , Δρ = 2000 kg m Excess mass, M  mass, M S  S 

= 100 m = 109 kg

 Note that:: at x = 0 m. ● gz has it’s maximum value directly above the sphere at x GM S  ● Th  The maximum value is gzmax =  z 2 of x where  where gz = (gzmax)/2 is called the half-width of the curve (x ½). ● The value of x Can show show that that x½ = 0.766 0.766 z ● Far away from the sphere, g z becomes very small

Gravity measurements are rarely made on a single profile. Usually they are made on a grid of points. This allows us to make a map of gz. Question: What will the map look like for the buried sphere

B2.2 Buried horizontal cylinder

When gravity measurements are made across a buried cylinder, it can be shown that the variation in gz will be.

 g  z  =

2Gπ a 2 z Δρ  ( x 2 + z 2 )

This curve is drawn below for a cylinder with radius, a = 50 m density contrast , Δρ = 2000 kg m -3

depth of axis, z axis, z = 100 m horizontal location, x location, x = 0 m

 Note that : of  g z  is located directly above the axis of the cylinder. ● the maximum value of g gzmax

=

2Gπ a 2 Δρ 

 z  ● From  From the the plots, plots, we can see that this value value is larger larger than than g zmax for a sphere sphere?? Why? Why? ● For a cylinder, we can show that the half-width x half-width  x½ = z 

Question: Compare the profiles across the sphere and a cylinder. Would this information allow you to decide if the buried object was a sphere or a cylinder?

Question: If gz  is measured on a grid of points, what will the resulting map look like? Would this be a better way to distinguish between a sphere and cylinder?

Forward and inverse problems in geophysics

B2.1 and B2.2 illustrate the gravity anomaly that we would expect to observe above a known geological target. This is called a forward problem  in geophysics, and is a useful exercise in understanding if measurements would be able to detect a particular structure. Forward problem:  Density model of Earth Earth > Predicted gravity data(anoma data(anomaly) ly)

However, we are usually more interested in solving the opposite problem. When gravity data has been collected in a field survey, we want to find out the depth and size of the target. This is called an inverse problem  in geophysics. Inverse problem :

 Measured gravity gravity data data > Density Density model of Earth

Example : Gravity data interpretation example Consider some gravity data collected on a profile crossing a spherical iron ore body.

● Where is the centre of the ore body?

 x

= _________metres

 g  z max = _________mgal

● What is the maximum value of g z?

● At what distance ( x ( x½) has gz fallen fallen to half this value? value?

x½

● The depth of the sphere can be derived using the equation

Rearranging this gives

= ________ _________me _metres tres

x½ = 0.766 z

z = 1.306 x½ z = ______metres

● To determine the excess mass, we use the equation max

We know z  and  and have measured  g  z

max

gz

=

GM S   z 2

 so we need to rearrange this equation to find  M s

 g  z max z 2

This gives

 M  s

Remember to convert g zmax from milligals to ms -2!

Ms = ________ _________kg _kg

=

G

● Can the radius ( a) and density contrast ( Δρ) of the sphere be determined?

B2.3 Uniform layer

A layer has a infinite extent, thickness ∆z and a density ρ. The gravitational attraction of this slab at P is:  g  z  = 2π Gρ Δ z   Note that that gz does not depend on the distance from the layer to the point P. Why? Consider the two density models shown below.

What can we say about the gravitational acceleration (g z) of the two models? This is an example of non-uniqueness  in geophysics, and occurs when more than one Earth model can explain the same set of geophysical data.

B2.4 General polygon – simple computer modelling

Simple density models such as a sphere and a cylinder can help us get an idea of what gravity anomaly (g z) is produced by simple density models. However, to represent the real Earth, more complicated model geometries are clearly needed. The development of  powerful, portable computers computers has has made this this possible. possible. In the lectures and labs we will use a very simple MATLAB program that computes g z for a N-sided prism. This is called 2-D modelling, since the polygon is assumed to be infinite in length in the strike direction (out of the plane of the page). Fully 3-D gravity modelling is also widely used, but will not be discussed in detail in Geophysics 224. B2.4 Example 1 – verification

When using a computer modelling program, it is important to be suspicious  about the answer. The first time the program is used, you should try and use a simple model for which the answer is already known. If you can get the same answer by two methods, then maybe you can trust the program. Consider using the program to compute the gravity anomaly of a cylinder. The cylinder can be approximated by an 8-sided polygon, as shown below. The dots show the exact (known) answer for a cylinder.

Question : Why is gz for the polygon less than for a cylinder?

With a 16-sided polygon, the agreement is significantly better.

B2.4 Example 2 – sedimentary s edimentary basin

An over simplified sedimentary “basin” is shown below sheet. Note that the “basin” is less dense than the basement rocks, so the pull of gravity is slightly weaker over the  basin. The horizontal line shows the expected value of gz  for an infinite layer with the same thickness and density as the basin.

 Now consider consider a similar similar basin, that is twice twice as deep deep and with with a lower density contrast. contrast.

 Note the following following two features: features: ● The minimum value of g z is the same as for the first basin. Why? ● The edge effects extend further into the basin than in the first case.

In the lab we will use this program to analyse some real gravity data collected across a sedimentary basin in Washington and an impact crater in Mexico.  MJU 2006 

224B3 Other factors that cause changes in g and need to be corrected  Note that gravity exploration is different to seismic exploration in the the following way: ●  In

a seismic survey, the travel time depends on just the velocity of the material on a path  between source and receiver.

● 

A gravity measurement is influenced by the gravitational attraction of local geological structure plus the mass distribution of the whole Earth, the moon, Sun, planets ….

B3.1 The Effect of latitude ●

The acce acceler lerati ation on of gravit gravity y at the the Equato Equator, r, g  E = 978,03 978,033 3 mgal mgal and and at the the poles poles g P = 983,21 983,219 9

mgal (Hammer, Geophysics, Geophysics, 8, 57, 1943). This difference is 5186 mgal, which is much  bigger than the anomalies we have considered, and needs to be accounted for in field measurements.

Three factors cause g to vary with latitude (A)The Earth is distorted by rotation

R E = 6378 6378 km km and and R P =6357 =6357 km. km. The The ratio ratio of flat flatteni tening ng is is appro approxim ximate ately ly 1/298 1/298.. Since Since a  point on the Equator is further  from the centre of the Earth than the poles, gravity will be slightly weaker at the Equator. We previously showed that for a sphere g ( r ) =

GM  E  2

r 

where the mass of the Earth, M E = 5.957 1024 kg. At the North North Pole, Pole, r = 6357 6357 km km and and g  P = 983, 983,219 219 mgal. mgal. If we move up 21 km to the radius of the equator, the decrease in gravity will be 6467 mgal Thus Thus g E = g P  - 6467 6467 mgal, mgal, which which is too too much much to explai explain n the the observ observed ed differ differenc encee betw between een the Equator and the Poles.

(B) Centrifugal forces vary with latitude The rotation of the Earth also causes gravity to vary. Imagine you are standing at the North Pole. The rotation of the Earth will not change  g, all that will happen is that you rotate once a day.  Now imagine you are at the Equator. If we could increase the rotation rate of the Earth enough, you would be ultimately be thrown into space (i.e. become weightless). Thus rotation makes gravity weaker at the equator. With With the the 1 rota rotatio tion n per per day day,, and and RE = 637 6378 8 km, km, can can show show tha thatt g  P = g E  +337 +3370 0 mga mgall Question : How fast would the Earth need to rotate to throw you into space?

(C)Mass distribution of the Earth The change in shape from a sphere to an ellipsoid redistributes the Earth’s mass. Thus results in more mass between points on the Equator and the centre of the Earth, than between the  poles and the centre of the Earth. This effect will make g  E > g P and is analogous the Bouguer correc correctio tion n we will will discu discuss ss in sectio section n B3.2. B3.2. Calcul Calculati ations ons show show that that g E ~ g P  + 4800 4800 mgal mgal

Combining these three effects (A,B and C) gives g P = g E + 6467 6467+ + 3370 3370 - 4800 4800 mgal mgal = g E  +503 +5037 7 mgal mgal ( appr approx oxim imat atel ely y as obse observ rved ed))

Equation for variation of g with latitude The variation of g from the Equator to the pole can be written as g (θ  (θ ) = 9.78031846 (1+ 0.0053024 sin ² θ  –  – 0.0000058 sin² 2θ  2 θ ) where θ  is  is the latitude in degrees. This equation is called the Geodetic Reference System  for 1967. More recent revisions are essentially the same, but with more decimal places.

Example : In Edmonton

g

θ =

53  30’ 25” N and the GRS67 equations gives ˚

= 9.7803 9.780318 1846 46 ( 1+0.00 1+0.00341 341790 7902-0 2-0.00 .00000 000539 5395) 5) m s = 9.81369388 m s

2

−

2

−

The variation of g with latitude is important when a survey extends over a significant northsouth distance. Differentiating the GRS67 equation with respect to θ  yields  yields

dg  d θ 

=

9.78031846 (0.0053024 x (0.0053024  x 2  2 sin θ cos θ - 0.0000058 x 0.0000058  x 4sin 2θ cos 2 θ)

=

0.04 0.0495 9526 26

ms

−

=

0.00 0.0008 0865 655 5

ms

−

=

86.550

mgal per degree

=

0.7868 0.7868

mgal mgal km

2

per per radi radian an

2

per per degr degree ee

1

−

( 1 degree degree latitu latitude de = 111 km) km)

All the these equations define the expected value of theoretical gravity (or normal gravity) at latitude θ. Differences between this value and what is actually measured are anomalies that we will analyse for information about subsurface density structure.

224 B3.2 How does elevation effect gravity measurements? (i) The Free air correction  Newton’s Theory of Gravitation states that at a distance, r , from the centre of the Earth g(r)

=

GM  2

r 

This means that as you move away from the centre of the Earth, the acceleration of gravity (g) decrea decreases. ses. In Edmo Edmonto nton, n, g = 9.81 9.81 ms gravity will decrease by Δg

= 3.086 Δh  x 10 = 0.3086 Δh

−

2

−

6

and if if you you move move up a distance, Δ distance,  Δh h, the acceleration of

2

−

ms mgal

Consider the exciting topography of a flat topped mountain:

Gravity measurements are made at points A and B. The difference in elevation means that g B will  be less than gB by  by an amount B

B

Δg

= 0.3086 x 100

= 30.86 mgal

When collecting gravity data, our real interest is to determine the density of the rocks below ground. The change in elevation from ‘A’ to ‘B’ will thus contaminate the data. The Free Air correction  is a mathematical way of undoing the effect of elevation. It allows us to correct the data collected at ‘B’ in order to make it equivalent to data collected at the same elevation as ‘A’. In gravity surveys, we always define a reference level  for the survey. Free Air corrections are made relative to this level. In general, any reference level could be chosen, but sea level is commonly chosen in coastal areas. In Alberta, the average level of the prairies would be a good choice. If a gravity measurement was made  Δh  Δh above the reference level, we must add C  FA = 0.3086 Δ 0.3086  Δh h

mgal 

C  FA is called the Free Air correction for a given gravity measurement. Similarly, if a gravity measurement was made  Δh  Δh below the reference level, we must subtract C  FA = 0.3086 Δ 0.3086  Δh h

mgal 

accurately must we know the elevation? Question : to keep data accurate to 0.1 mgal, how accurately

(ii) The Bouguer correction

Unfortunately, this is not the end of story! Compare the gravity measurements at ‘A’ and ‘B’. At  point A, the gravity measurement is solely due to structure below the reference level (blue). At ‘B’ the gravity measurement is due to structure below the reference level, plus the gravitational  pull of the 100 metres of mountain (red). (red). The net result is that g B > gA B

From section B2.3 the magnitude of this extra gravitational attraction is approximately gB - gA = 2πG ρ  Δh  Δ h where ρ where  ρ is  is the density of the mountain. Thus to remove this effect we need to subtract

C  B = 2πG

from from the the obser observe ved d grav gravity ity

measurement at ‘B’. This is called the Bouguer correction and C B = 0.00004193  ρ Δ  ρ  Δh h  mgal  Note that to apply the Bouguer correction we need to estimate  ρ,  ρ, the density that lies between ‘B’ and the reference level. Using the value  ρ =  ρ = 267 2670 0 kg kg m

3

−

this this give givess

C  B = -0.1119  Δh  Δh  mgal This value represents an average density for crustal rocks. Other information ( e.g  borehole  borehole gravity data or Nettleton’s method) method) may be used to give a better estimate of the density.

Pierre Bouguer

Summary Measurement above reference level

Add Free Air correction

Subtract Bouguer correction

Measurement below reference level

Subtract Free Air correction

Add Bouguer correction

MJU 2006

224 B4. Instrumentation and field procedures

224B.4.1 Instrumentation for gravity exploration - gravimeters 4.1.1 Absolute measurements of g

A range of devices can be used including reversible pendulum and free fall devices. They are generally more expensive and slower to operate than a mass-on-a-spring gravimeter. For details of instrumentation Application to tectonics in BC

http://www.microgsolutions.com/ http://www.pgc.nrcan.gc.ca/geodyn/abs_grav.htm http://www.pgc.nrcan.gc.ca/geodyn/cascadia.htm

4.1.2 Relative measurements of g (measure differences in g)

Since we have seen that it is the differences in gravity that is generally more important than absolute values, we do not need absolute measurements of gravity at every survey location. Often relative gravity measurements can be made over a survey area, and then tied to an absolute value by using the relative gravimeter at a location that was  previously surveyed with an absolute gravimeter. (a)  Portable pendulum Accuracy around 0.25 mgal when popular in the 1930’s (b) Mass-on-a-spring (b)  Mass-on-a-spring gravimeters The mass experiences a force of  F = mg  and the spring stretches an amount  s.  s. Hooke’s Law states that  F = ks = mg where k  is   is a measure of the stiffness of the spring (the spring constant). If the gravimeter is then taken to a location where the acceleration of gravity is stronger by an amount δg, then the spring will stretch a little bit more, δs. k(s+δ s) = m(g+δ g) Subtracting these two equations gives δ g

δs

=

k δ   δ s    m

can be very small, so various engineering features are used to amplify the movement. Typical accuracy is 0.01 mgal. mgal. What change in elevation produces 0.01 mgal?

Two widely used gravimeters are: ●LaCoste-Romberg ●Worden

gravimeter ( http://www.lacosteromberg.com http://www.lacosteromberg.com))

gravimeter (http://www.gravityservices.com/meter.htm (http://www.gravityservices.com/meter.htm))

4.1.3 Gravity gradient measurements (dg/ dx dx)

These were initially made with torsion balances. In a torsion balance, two masses are suspended from a beam. The force required to twist the thread is determined by calibration. The horizontal pull of gravity will tend to orient the beam so it points towards or away from concentrations /deficits of mass. The Cavendish torsion  balance allowed measurements of G  to be made. An alternative design was developed by Baron von Eötvös in Hungary in the later 19 th century and was used in the first geophysical discovery of oil in the early years of the 20 th century. The torsion balance was quite effective at locating salt domes, owing to the lower density of the salt compared to the surrounding sedimentary rocks. The torsion  balance required a light-proof and wind proof chamber to house house the instrument.

Figures from Applied from Applied Geophysics, Geophysics, Eve and Keys, Cambridge University Press, 1928

4.1.4 Marine gravity surveys

Conventional gravity measurements can be made by stopping the ship and lowering a gravimeter to the sea floor for each measurement. However, this is very time consuming, especially in deep water.

To measure gravity on a moving ship, the gravimeter is placed on a gyroscopically stabilized platform and for each measurement, g is measured for several minutes. The acceleration due to the waves is equally positive and negative and thus the average value is due to the acceleration of the Earth. This steady velocity  of the ship must also be accounted for through the Eötvös correction. Typical accuracy of marine gravity data is around 1 mgal. 4.1.5 Airborne gravity surveys

The problems encountered in marine gravity surveys are even more severe in airborne exploration (turbulence can produce an acceleration greater than g)  and velocities are higher. However, airborne measurements are excellent at covering large, roadless areas very easily. In contrast to marine surveys, the acceleration of the aircraft can be independently measured through accurate measurements of aircraft elevation as a function of time. The acceleration of aircraft can be computed and then removed from the acceleration recorded on board.

From, Sigmund Hammer, Airborne gravity is here, Geophysics, Geophysics, 1983. Note that: ● Upward continuation is a mathematical method that uses gravity data at one elevation to compute the gravity data that would be measured at a higher elevation. ●  Gravity measured at 300 m elevation in a helicopter agrees with the upward continuation of the surface gravity measurements (land Bouguer). ●By moving further away from the target, small features such as the (higher density) cap rock are not detected. This is a consequence of the inverse square law.

224 B.4.2 B.4.2 Instrument drift and tidal effects Instrument drift If repeated gravity measurements are made at the same location  with a relative gravimeter, the answer will change over time.

Occasionally, this is due to changes in the density of the Earth. However, it is usually due to the spring in the gravimeter stretching over time. This effect can be removed by establishing a base station that is visited repeatedly during the survey. This problem arises because the gravimeter cannot be in two places at once. For example, when the first measurement is made at Base 1, we would like to know the gravity measurement that would be made at the same time at IGSN71 base. This cannot be done, so linear time variation is assumed between the two measurements at the IGSN71 base. This allows us to calculate the difference in gravity measurements (dial reading) between Base 1 and IGSN71 base. Note that the Dial reading is a relative measure of g (mgal).

(1) IGSN71 g = 981,000 mgal (2) Difference between IGSN71 and Base 1 is 9 mgals (3) Thus, at Base 1, g 981,009 mgals

(absolute) (relative) (absolute)

Tidal effects In additional to instrument drift, gravity measurements made at the same location will vary with time due to tidal effects. This is dominantly due to the gravitational pull of the sun and moon. The magnitude can generally be predicted accurately and removed from the measured gravity data

224B.4.3 Gravity survey procedures gravity data on a 2-D grid and repeat measurements at cross-over points where lines intersect. This will give a good idea of the repeatability of the measurements.

● Collect

● Set

up a series of base stations to remove instrument drift and tidal effects from data. These should be visited several times a day.

can obtain absolute  gravity relative measurements from a gravimeter by making measurements at pre-surveyed stations  where the absolute values is already known.

●We

A network of these stations has been established across Canada and is called the Canadian Standardized Gravity Network (CGSN).

Details can be found at http://www.geod.nrcan.gc.ca/index_e/products_e/passiveNetworks_e/cgsn_e/cgsn_e.html ● Details

of a base station at the International Airport are listed on the following page and can be found at: http://www.geod.nrcan.gc.ca/database7/grav/sd/stage/95001993EN.html

Here the value is 981117.890 981117.890 mgal Next time you have time to spare waiting for a flight, see if you can find the benchmark.

Canadian Gravity Standardization Network (CGSN) STATION IDENTIFICATION (IGNS71) EDMONTON ALB 9500-1993 Gravity: 981117.890 mgal Last Inspection: 08/1996

STATION COORDINATES (SCALED) Latitude :  N 53° 18' 29" Longitude : W 113° 35' 1" Elevation : 723 m

STATION INFORMATION AND LOCATION

The station is located in the Terminal of the Edmonton International Airport at Nisku. It is on the floor, in the corner, E of the Information counter on the ground floor of the terminal  building. The station is monumented monumented by an aluminum disc. SKETCH:

DISTANT VIEW:

CLOSE-UP VIEW:

Geodetic Survey Division, Ottawa http://www.geod.nrcan.gc.ca/index_e/products_e/passiveNetworks_e/cgsn_e/cgsn_e.html

●Often

several survey crews needed for each gravimeter crew. It is vital to know the elevation of each measurement location. Differential GPS is sometimes good enough,  but levelling may be needed.

●  Marine

and airborne surveys. Ground is covered much more quickly than with land based methods, but in general observations made further away from potential targets.

●Satellite

gravity: Cannot actually measure acceleration while the satellite is in free fall. However, from detailed observations of the orbit, or by measuring the altitude of the satellite above the Earth’s surface, subsurface density structure can be inferred.

Sources of error in gravity data

(a) Vertical position

0.3 mgal per m

(b) Horizontal position

0.0008 mgal per m north-south

(c) Inadequate or incorrect terrain corrections

~ 1 mgal or more

(d) Instrument drift

0.5 mgal per month

(e) Incorrect density used in Bouguer correction e.g. change in elevation = 100 m, density error 0.1 g cm

 MJU

Geophysics 224 January 2006

−3

gives an error of 0.42 mgal

Geophysics 224 B5 : Interpretation of gravity data B5.1 Bouguer anomaly maps ●

  Once the Free Air and Bouguer corrections have been made, the Bouguer anomaly should contain information about the subsurface density alone. The effect of latitude and elevation should have been removed. A map of the Bouguer anomaly gives a good impression of subsurface density. Low (negative) values  of Bouguer anomaly indicate lower density  beneath the measurement point. High (positive) values   of Bouguer anomaly indicate higher density beneath the measurement point.

From M. Pilkington et al, Can. J. Earth Sciences, 37, 1453-1471, 2000. Note decrease in Bouguer gravity anomaly towards the southwest. This is due to (a) thickening of the Alberta Basin and (b) thicker crust beneath Rocky Mountains.

B.5.2 Regional trends and residuals  Measurements of gravity at the Earth’s surface depend on Earth structures structures ranging from spatial scales of 1m to 10,000 km.

●

●

Deeper density structures produce gravity anomalies with a long spatial wavelength. For a body such as a sphere, this would be expressed as a large value for the half-width (x ½) .

●

  Similarly shallow density structures wavelengths or small half-widths.

produce gravity anomalies with short spatial

●

  Often the longer wavelength effects are called regional trends, while the shorter wavelength features are called residuals or  anomalies   anomalies. The distinction between the two is somewhat arbitrary.

●

  In shallow gravity exploration, we are generally interested in the short wavelength gravity anomalies. The long wavelength regional trends can make it difficult to analyse the short wavelength residuals. Thus we need to find a way to remove the regional trends and emphasize the anomaly more clearly. In other studies the regional trend could be what we wish to study (one person’s noise is another persons signal!)

●

 To emphasize shallow structure, the regional trend can be computed and then subtracted from the measured Bouguer anomaly data. However care is needed not to eliminate useful data during this process.

●

  Sometimes the regional trend can be approximated by a straight line. In other case it should be approximated with a more complex mathematical form.

density structures (a) higher density bedrock that dips to to Example 1 : This shows two density the right and which produces a slow decrease in g B from left to right. (b) a shallow  buried cylinder. It is difficult to determine the half-width of the anomaly due to the cylinder in the presence of the linear regional trend.

●

The dotted line denotes the regional trend, and is obtained by finding the best-fitting straight line to the data. When the regional trend is subtracted from g B the anomaly of the cylinder is much easier to interpret. B

Example 2: This features the same dipping bedrock layer and a shallow, low density river channel. Again, removal of the regional trend makes the effect of the river channel much clearer. ●

●

Example 3: Real gravity data collected in a geothermal exploration project in the Philippines (courtesy of PNOC). These gravity data were analysed in Lab 1 after the regional trend was removed.

●

 The above examples are for a single profile. However, regional trends are often removed from Bouguer anomaly data defined on a grid (i.e. a Bouguer anomaly map).

●

 Must be careful not to remove useful data when removing regional trends from the data.

B.5.3 Half width techniques  For a number of shapes the width of the anomaly can give an approximate depth to the density anomaly. e.g. sphere : depth of centre = half-width cylinder : depth of centre = 1.304  x half-width

●

 However in some cases this approach can be very wrong! Dike model with MATLAB program, Anticline ( Geophysics, 12, 43, 1947)

●

B.5.4 Two-dimensional computer modelling  This allows the user to handle a more realistic density model. The data is usually fit with a trial and error approach. This is illustrated with the simple MATLAB program used in Laboratory 2 to model a couple of real datasets.

●

B.5.5 Horizontal derivatives of gravity Taking the horizontal derivatives of the Bouguer anomaly emphasizes changes in the horizontal gradient. This is effectively another way of removing (or suppressing) the regional trends in the data.

●

 Consider the application to the basin example shown below. The gradient gradient is very effective at -9 -2 defining the edge of the basin. The unit E is the Eotvos. 1 E = 10 s

●

 In the example below, note that the second derivative has a different sign depending on the geometry of the edge of the basin..

●

Gravity gradients : Alberta basement

The primary domains of basement geology in Alberta are shown on the left. Note that the several domain boundaries are coincident with high values of the horizontal gravity gradient (e.g. Vulcan structure in Southern Alberta and Snowbird Tectonic Zone). These variations are clearer than in the Bouguer anomaly map shown earlier in this handout.

Gravity gradients : Chicxulub Impact Crater, Mexico

The horizontal gradients of the Bouguer gravity anomaly data from the Chicxulub impact structure defines the edges of the structure very clearly. Data from Alan Hildebrand, University of Calgary. Noise in field data can make this procedure unstable, so smoothing may be needed.

Geophysics 224 MJU 2006

Geophysics 224 B6 Applications of gravity exploration B6.1 Location of caves and man-made underground cavities

●Cave

location in karst terrain. Taken from Burger Figure 6-37. The caves  produce a decrease in the Bouguer gravity anomaly. ●  Sand

and clay have a lower density than limestone. Note that the variable thickness of sand and clay can mask the effect of the voids. ● The

-30 mgal anomaly looks very big for such a small cave! Can you verify verify the result with the formula for a cylinder?

●  Microgravity

was recently used by United Nations weapons inspectors to look for underground bunkers in certain Middle Eastern countries.

“Microgravity meters -- also called gravimeters -- measure minute differences in gravitational pull at one site versus another. Large underground voids, such as tunnels or weapons  production facilities, slightly lower Earth's gravitational pull at the surface right above the voids. Gravimeters can detect these differences, indicating where such facilities might exist. According to a source familiar with the inspections, gravimeters operate too slowly to efficiently scan large areas. However, they work well within a single structure, such as a palace or a bunker, where single and/or multiple basements are suspected” (Geotimes, 2002)

Figures from an interesting article by Won et al., (2004). See also http://www.agiweb.org/geotimes/nov02/WebExtra112702.html ● Kearey

Figure 6.30 Microgravity detection of cavities at the site of a proposed cooling tower for a power station. Note that contours are in gravity units (g.u.).

B6.2 Geotechnical studies ● Size

and geometry of landfills. The contents of the landfill are usually lower in density than the rock that was removed. Burger Figures 6-38 and 6-39 show a study of a known landfill that was used to see what gravity anomaly was associated with a small landfill. Before and after air photos were used to estimate the depth of the landfill.

● What ● How

is the average density of this landfill? accurate must the vertical surveying be for data that is accurate to 0.1 mgal?

● Other

applications of gravity exploration in hydrogeology. Kearey Figures 6.28 and 6.29 show how gravity data were used to detect buried river channels containing aquifers near Antofagasta in Northern Chile. The gravity data do not detect the aquifers directly, but locate structures that contain the aquifers.

● Depth

to bedrock studies

B6.3 Mineral exploration

Ore bodies are often higher density than the host rock and can produce positive gravity Bouguer anomalies. Gravity can be used to estimate the excess mass of an ore deposit, using Gauss’s theorem. While While non-uniqueness prevents the spatial distribution being uniquely determined, the total excess mass  can be estimated reliably (see Kearey 6.10.3).

The example below comes from Voisey’s Bay, Labrador. This massive sulphide deposit has a pronounced pronounced positive gravity gravity anomaly. Note that gravity inversion is an automated procedure that determines a density model that fits the measured gravity data. An inversion is a solution of the inverse problem, and non-uniqueness must be taken into account.

More information : http://www.geop.ubc http://www.geop.ubc.ca/ubcgif/caseh .ca/ubcgif/casehist/voisey/intro.html ist/voisey/intro.html

B6.4 Hydrocarbon exploration

Torsion balances were used in the 1920’s in the first discovery of an oil deposit in the United States with a geophysical method (Nash Dome, Louisiana). The torsion  balance was quite slow and was rapidly replaced by pendulums in the 1930’s. Applications of gravity exploration in hydrocarbon exploration include: ● Reconnaissance

surveys to map the thickness and extent of sedimentary basins. Can also locate buried rift structures and basement faults. Generally basins are gravity lows. Reef trends can sometimes be identified.

In addition to surface and airborne gravity data, satellite gravity data  has contributed to exploration in offshore areas.

Figures from: http://topex.ucsd.edu/marine_grav/mar_grav.html ●

Salt structures. Salt has a lower density than surrounding sedimentary rocks, and gravity data can be used to characterize the geometry of salt layers and domes. Locating the base of a salt body  is difficult with seismic reflection data, and alternative geophysical methods are being used. In combination with seismic reflection data, gravity can be used to give a joint interpretation. This can include hypothesis testing of a velocity model. Since gravity data is usually available on a grid, it can be used to extrapolate features between, or beyond, seismic coverage.

●Continental

scale data gravity databases can also be used. More effective when regional features due to isostatic effects are removed (Chapin, 1998 in volume edited by Gibson and Millegan, 1998).

B6.5 Regional scale studies of lithospheric structure Impact craters : The Chicxulub Impact crater in Mexico was initially discovered from Bouguer anomaly data. The ringed shaped gravity low is caused by a basin that is filled with lower density sedimentary rocks that were deposited in the crater after the impact. The Bouguer gravity high in the centre is due to crystalline  basement rocks that rebounded immediately after impact. In Lab 2 we will interpret these gravity data and estimate the thickness of the basin.

●

●

Isostacy : Gravity measurements have also led to the concept of isostacy. Based on observations by Pierre Bouguer in South America (1735), and George Airy and John Pratt in India (1800’s) who noted that pendulums were not deflected as much as expected by large mountain ranges. These observations can be interpreted to show that these mountains and plateau are supported by a low density root. This is very similar to the way that an iceberg floats in seawater, where most of the ice is below sea level.

●

Basin structure and earthquake hazards : The Seattle Fault forms the Southern edge of the Seattle Basin, and may have caused M=7 earthquakes in the recent geological past. In Lab 2 we will model gravity data that have been used to study the geometry of this fault. In this study data were collected along the shore of Lake Washington. This allows measurements to be made at a constant elevation and removes the need for Bouguer or Free Air corrections to be made.

Modelling of these gravity data also provides a useful illustration of the inherent non-uniqueness associated the interpretation of gravity data.

B6.6 Temporal variations in gravity

Time variations in gravity have been observed on active volcanoes as magma intrudes prior to an eruption. What could cause an increase or decrease in gravity in such a case?

The data above were collected at Long Valley Caldera, California. It is not clear if this will provide a useful tool for predicting eruptions. Typical changes are of the order of 50 μgal. Data analysis requires that time variations in elevation must also be accounted for. Why is this important?

References

M. Battaglia et al, Magma intrusion beneath Long Valley Caldera confirmed by temporal changes in gravity, Science, 285, 2119-2122, 1999. Geologic Applications of gravity and magnetics: Case Histories, edited by R.I. Gibson and P.S. Millegan, Society of Exploration Geophysicists, 1998. Gravity and magnetic methods at the turn of the millennium, Geophysics, 66, 3637, 2001. Won, I., et al., Geophysics and weapons inspection, The Leading Edge, 658-662, July, 2004

 MJU 2006 Geophysics 224 224

Geophysics 224 - 2006 B. GRAVITY EXPLORATION - Outline 1. Introduction 1.1  Newtonian gravitation for a point mass 1.2

Density of rocks and minerals

1.3

Units for gravity measurements (milligals)

1.4

Approximate calculation to estimate size of gravity anomalies

2. Gravity anomalies of some simple shapes 2.1 Buried sphere 2.2 Horizontal cylinder 2.3 Uniform layer  2.4 General polygon – simple computer modelling

3. Other factors that cause changes in g 3.1 Latitude of measurement location ● The Earth is rotating and not exactly spherical. ● GRS67 formula 3.2 Elevation of measurement location ●Free air and Bouguer corrections

4. Instrumentation and field procedures 4.1 Measuring devices (gravimeters) ● Absolute gravity measurements : free fall devices and pendulums ● Gravity gradient measurements ● Relative measurements : pendulums, mass-on-a-spring gravimeters ● Airborne and marine measurements now possible. 4.2 Survey procedures ● Set up base stations. Tie to CGSN stations to give absolute values. 4.3 A few more corrections! ● tidal effects, instrument drift

5. Interpretation of gravity data 5.1 Bouguer anomaly maps 5.2 Regional and residuals 5.3 Half width techniques 5.4 Forward modelling with a general 2-D model 5.5 Use of horizontal gravity gradient to define edges

6. Applications 6.1 Locating caves and voids. Geotechnical, archaeological and military applications. 6.2 Geotechnical and hydrogeology applications. Definition of landfill geometry. 6.3 Mineral exploration. Ore bodies are usually high density. 6.4 Hydrocarbon exploration. Eötvös torsion balance used in early 20 th century to locate salt domes. Airborne and marine data used in reconnaissance of sedimentary basins. Gradient techniques for imaging base of salt layer. 6.5 Regional scale studies of crustal and upper mantle structure. e.g. Western Canada, Chicxulub impact crater in Mexico. 6.6 Temporal variations in gravity

 MJU 2006 

224 C1 Basic principles of seismology C1.1 General introduction to wave phenomena A wave can be defined as a periodic disturbance  that transmits energy through a medium, without the permanent displacement  of the medium. Also required that energy is converted  back and forward between two different types. Before considering the details, it is useful to review some basic definitions. Consider the two waves shown in the MATLAB movie waves.m

 Frequency (f)  (f)  : The number of cycles a given point moves through in 1 second. Frequency is measured in Hertz (Hz). If the frequency is very low, then it is common to refer to the period (T) of the signal in seconds. T = 1/f  Angular frequency ( ω ): Frequency is the number of rotations per second. The angular frequency is the number of radians per second and given by ω=2π  =2π  f  Wavelength (  λ ) :  ) : Distance in metres between two points of the wave having the same phase ( e.g. two crests or two troughs). If the waves moves at a velocity, v, then v = f λ f  λ ● Can you check that this relationship is correct for the figure above? ● Note that points on the wave move up and down, they are not translated to the right. ● In this case, the velocity is independent of frequency. This type of wave behaviour is called non-dispersive. If velocity varies with frequency, the wave is said to be dispersive

1

In seismology, we need to understand how waves will travel in the Earth. For example, how fast will they go, which direction, how will amplitude vary with distance etc. In general this requires the solution of some complicated differential equations. In Geophysics 224 we will approach this subject through visualization. Wave propagation can be considered in two ways, by considering either wavefronts or rays. These are complementary ways of talking about waves:

Rays  denote the direction in which the wave travels. Wavefronts  are points on the wave with the same phase (e.g. a line along the crest of a wave is a wavefront).

 Note that wavefronts and rays are at right angles to each other.

One way to visualize wave propagation over time is through Huyghens Principle. This states that all points on a wavefront can be considered secondary sources of wavelets. These secondary wavelets propagate outwards and at a time later, the overall wavefront is the envelope of secondary wavelets. Examples for a point point source is shown below.

2

C1.2 Wave propagation in the Earth Having considered some general aspects of wave propagation we now need to consider how waves will propagate in Earth. Waves in the Earth can be divided into two main categories: (a) Body waves travel through the  bulk medium. (b) Surface waves are confined to interfaces, primarily the Earth-Air interface.

C1.2.1 Body waves

Body waves in the Earth can be divided into two types:

P-waves : Particle motion is in the same direction as the wave propagation. They are also called compressional or  longitudinal  longitudinal waves and

S-waves: Particle motion is at right angles to the wave propagation. Also called shear waves or transverse waves) The velocity of body waves can be calculated from the properties of the material, as outlined  below. More M ore details are given by Kearey pages 22-23, 3 rd Edition. The mathematics is simplest for the geometry shown shown below for a column of rock rock with cross-sectional area ( A).  A).

3

If a force ( F ) is applied at the left end, this can be considered an applied longitudinal stress =  F/A.  F/A. The rock deforms elastically in response to the applied stress. The deformation of the leftmost disk can be quantified as the longitudinal strain = strain =

Δ L  L

The behaviour of the medium is described by the longitudinal modulus,

ψ = ψ  = longitudinal stress / longitudinal strain

which is a measure the stiffness of the rock (i.e. ( i.e. how  how much stress is needed to produce a given strain). This strain produces a force that will cause the the shaded section of the rock to accelerate to the right. This lowers the stress to the left, but increases it to the right. This causes the next section of the rock to move and so on. Can show that a wave motion will move down the column at a velocity

v=

ψ   ρ 

where  ρ is  ρ is the density of the material. Note that the stiffer the medium (larger ψ ) the greater the force on the shaded cylinder, thus acceleration is higher and wave velocity is greater. Similarly, as density increases, the shaded section becomes heavier and it’s acceleration (and wave velocity) for a given force will decrease. In general, the calculation of velocity is more complicated as the deformation will involve both compression and shearing. Two other modulii must be defined to fully fully understand these effects.

Bulk modulus, K modulus, K = volume stress, P volume strain, Δ strain,  ΔV/V  V/V  and Shear modulus, μ = shear  = shear stress, τ   shear strain, tan θ 

4

It can be shown, that these modulii are related as ψ  =  K  +

4 3

μ 

Thus the P-wave velocity can be written as 1

⎡ K  + 4 μ ⎤ 2 3 ⎥ v P  = ⎢ ⎢  ρ  ⎥ ⎣ ⎦ and the S-wave velocity as 1

⎡ μ ⎤ 2 v S  = ⎢ ⎥ ⎣ ρ ⎦ ● Note that P-waves always travel faster than S-waves ●  In a liquid  μ=0  μ=0 while  K   is always non-zero. Thus only P-waves can travel in a liquid, since shear stresses cannot exist. ● These expressions for v P  and vS  do not depend on frequency, thus body waves (both P-waves and S-waves) are non-dispersive. ● Note that as the rock cylinder is stretched, it will get longer and thinner. This effect can be quantified through Poissons ratio. This is defined as: σ = σ  =

lateral strain longitudinal strain

● Note that individual values of v P  and vS depend on several modulii and density. This can make it difficult to compare the velocities of similar rocks. It can be shown that

v p /v s = [2(1-σ  [2(1-σ  )/(1-2σ   )/(1-2σ  )]½

For typical consolidated crustal rocks, σ ~ 0.25 and v p /v s ~ 1.7. An increase in v p /v s and/or Poisson’s ratio can be indicative of the presence of fluids.

5

C1.2.2 Surface waves These waves are localized at the Earth’s surface and can be divided into two types. The particle motions are illustrated in Kearey Figure 3.4 Rayleigh Waves are like ripples on a lake. Particle motion is in a retrograde ellipse.

motion analogous to S-waves. Love waves have a horizontal particle motion

 Note that: ●The velocity of a Rayleigh wave does not vary with frequency when travelling in a uniform medium and it is slower than an S-wave. In a layered Earth the velocity of a Rayleigh wave varies with frequency (it is dispersive) and can be used to infer velocity variation with depth. ●Love waves only exist if the Earth is layered and are always dispersive. ●Large earthquakes can generate surface waves that travel around the globe. They can be large in amplitude and cause a lot of damage during earthquakes. ●In exploration seismology, ground roll  is a Rayleigh wave that travels across the geophone array.

6

C1.3 Typical seismic velocities for Earth materials Typical values for P-wave velocities in km s -1 include: Sand (dry) Wet sand Clay

0.2-1.0 1.5-2.0 1.0-2.5

Tertiary sandstone Cambrian quartzite Cretaceous chalk Carboniferous limestone Salt

2.0-2.5 5.5-6.0 2.0-2.5 5.0-5.5 4.5-5.0

Granite Gabbro Ultramafics

5.5-6.0 6.5-7.0 7.5-8.5

Air Water Ice Petroleum

0.3 1.4-1.5 3.4 1.3-1.4

● why does v p apparently increase with density? e.g. for the sequence granite-ultramafics. 1

⎡ K  + 4 μ ⎤ 2 3 ⎥ suggests that v p should decrease as density increases. The equation v P  = ⎢ ⎢  ρ  ⎥ ⎣

⎦

 P-wave velocity of porous rocks A rock sample comprises a matrix of rock grains with pores in between. The fractional volume occupied by the pore space is termed the porosity (Φ ( Φ). Pores are generally filled with a fluid such as air, water, oil or gas. If the matrix has a density  ρm  and the  pore fluid has a density  ρ f , then the overall density of the rock will be  ρ =  ρ = ρ  ρ f Φ + (1- Φ ) ρ  )  ρm Thus if the velocities in the matrix and fluid are vm and v f  it is quite simple to show that

1 v

=

Φ v f 

+

(1 − Φ) vm

This equation is sometime called the timeaverage equation. equation. An example showing how porosity controls v P   is in Kearey Figure 3.6. Note that the pore fluid is water, as evidenced by the density at 100%  porosity.

7

Example 1 : A rock has 30% porosity and the P-wave velocities in the rock matrix and pore fluid are 1.5 and 2.9 km s -1 respectively. What is the overall seismic velocity? Example 2 : Consider the gas reservoir shown below. The reservoir sand has 10% porosity. The following velocity data data was obtained from well-log information

 Material Sandstone Water Gas Shale

v P  (km s-1 ) 4.27 1.50 0.3 2.4

What is the average P-wave velocity in the water-saturated reservoir? What about in the gas saturated reservoir? Example 3 : In a similar reservoir, the P-wave velocity in the gas filled horizon was measured to  be 2.2 km s-1. What porosity does this imply?

Variation of seismic velocity with depth

With increasing depth, compaction increase the density of a rock through reduction of pore space. The rigidity of the rock also increases with depth. The net effect is that velocity will increase with depth, even if the lithology does not change.

The increase in velocity with burial depth is quite well defined, and deviations from the expected velocity at a given depth are sometimes used to infer uplift histories in basins. This technique is  based on the observation that the compaction is essentially an irreversible process (Telford Figure 4.21)

8

Rippability

In geotechnical studies, a knowledge of P-wave velocity can be used to determine if a  particular rock type can be ripped with a given tractor. This also gives a useful  Kearey Figure 5.25) compilation of the range of velocities for a given material. ( 

As a seismic wave travels through the Earth, several factors will change the direction  and amplitude of the waves. When detected at the surface, an understanding of these factors can tell us about sub-surface structure. In sections 1.4 – 1.7 we will examine a number of these factors.

9

C1.4 Reflection coefficients at normal incidence Consider a seismic wave that is travelling vertically downwards, as shown in the figure below. If the wave reaches an interface, some of the energy will be transmitted, and some will be reflected.  Note that the waves travels from medium 1 into medium 2. For each medium, the impedance is defined as the product of density ( ρ) and seismic velocity (v), Z = ρv

The reflection coefficient, R coefficient, R,, is defined as

 Ar   Ai

=  R =

v 2 ρ 2 − v1 ρ 1 v 2 ρ 2 + v1 ρ 1

=

 Z 2 − Z 1  Z 2 + Z 1

and the transmission coefficient, T , is defined as

 At   Ai

= T  =

2v1 ρ 1 v 2 ρ 2 + v1 ρ 1

=

2 Z 1  Z 2 + Z 1

The reflection coefficient is a measure of the change in impedance across the interface. These equations are called the Zoeppritz equations. If the wave is incident at an angle the become more complicated. Sometimes reflection co-efficients are expressed in terms of energy, not amplitude. The equations for this situation are listed in Kearey, chapter 3.

10

Example 1 The MATLAB script reflect_v1.m generates a movie showing how a seismic pulse propagates in the Earth. For simplicity, density is constant and only velocity varies with depth. In this example there is an increase in seismic velocity (and impedance) at 500 m. Note that the reflection has the same polarity as the down going pulse. What are the values of R of  R and  and T ? Is the largest amplitude in the reflected or transmitted wave?

Example 2 Same geometry as Example 1, but a decrease in velocity (impedance) occurs at 500 m depth.  Note that the reflected pulse has the opposite (negative) polarity t o the down going pulse. Is this what the Zoeppritz equations predict? What are the values of R of  R and  and T ? Is the largest amplitude in the reflected or transmitted wave?

11

Example 3 This has two interfaces, each associated with an increase in velocity. Both arrivals have a positive polarity , showing that they originate at interfaces where the velocity (impedance) increases. What are the reflection co-efficients at each interface? What other factor makes the amplitude of the second reflection smaller than the amplitude of the first?

12

Example 4 This has a low velocity layer (LVZ). The first reflections has negative polarity show a decrease in velocity, while the second has positive polarity, showing an increase in seismic velocity. Note that a pulse reverberates  within the low velocity layer producing multiple later arrivals, with ever decreasing amplitude.

Example 5 Consider a gas reservoir reservoir with 10% porosity where where the rock matrix has a v p = 3000 ms-1. The overlying rock has v p = 2800 ms -1. ●What is velocity in the reservoir? (use time average equation) equation) ●What value of R is expected for a reflection from the top of the reservoir? ●What value of R is expected where the 2800 ms

-1

-1

and 3000 ms  layers are in direct contact?

This is an example of a bright spot, which is a high amplitude reflection. In the Earth, typical reflection coefficients are ±0.2 with the maximum values ±0.5

13

C1.5 Reflection and refraction at non-normal incidence In a more general case, the seismic wave will be incident on an interface at some angle of incidence, θi. Note that the angle of incidence is measured from the normal to the ray. Snell’s Law was developed for optics, but can equally be applied to the seismic case. Consider a P-wave that strikes the interface shown below. In this case v 2 > v1

Reflected and refracted P-waves are generated from the incident P-wave. For the reflected Pwave, θ r r = θ i . Snell’s Law states that the ray parameter,  p,  p, for the incident and refracted waves   will be constant.

 p =

sin θ i v1

=

sin θ t  v2

=

sin θ r  v1

Thus with a velocity increase (v2 > v1) the ray is refracted away from the normal. If the velocity decreases (v  (v2 < v1) then the ray is refracted towards the normal.  Note that if v2 > v1 then there will be a value of θ i which results in sin in  sin θ t t = 1. This gives a value   of θ t t = 90  and the refracted waves travels horizontally. In this configuration θ i = θ c and is called   the critical angle. ˚

sin 90 o v2

=

sin θ c v1

and by rearranging we find that

⎛  v1  ⎞ ⎟⎟ v ⎝  2  ⎠

θ c = sin −1 ⎜⎜

The wave travelling horizontally is called a head wave. For a head wave to develop, we must have v2 > v1. Using Huyghen’s Principle, it can be shown that the head wave will generate upward propagating wave at an angle θ c to the normal. When these waves reach the geophones they are called refracted arrivals. When θ i > θ c the wave is totally reflected.

14

In general we must also consider that the incident P-wave will generate a reflected S-wave and a refracted S-wave through a process called mode conversion. Again Snell’s Law can be used to calculate the angles of reflection and refraction.

 p =

sin θ i v1

=

sin θ  P  t  v P  2

=

sin θ t S  v 2S 

The amplitudes of the transmitted and reflected waves can be calculated from Zoeppritz’s equations, which have a more complicated form than those listed in C1.4 (see  Applied Geophysics, Geophysics, Telford, 2 nd Edition, page 157 for details and some examples). The amplitude of the reflected wave varies with angle, and is sensitive to the change in impedance across the boundary. This is the physical basis of the amplitude versus offset (AVO) technique that we will discuss later on.

C1.6 Diffraction Seismic energy can sometimes travel in regions where ray theory (Snell’s Law) does not predict that it will go. This is typified by diffraction which occurs when a wave strikes an object that is significantly smaller than a wavelength. A diffractor effectively radiates seismic energy in all directions. When a wave strikes a corner, Huyghens Principle shows that the corner will generate waves that  propagate in all directions. We will see in real seismic reflection data that these waves can have significant amplitudes and are detected over a wide area at the surface. The mathematical process of migration  is used to trace the pattern of diffracted waves back the point source where they originated. This allows subsurface geometry to be correctly imaged in 2-D and 3-D.

15

C1.7 Factors that cause the amplitude to change as wave propagates (a) Geometrical spreading Imagine a wave travelling outward from a point source. If the wave has travelled a distance r  ,  , 2 then the wavefront covers an area A area  A =4π  =4π r  r  .At this point the wave has an amplitude of X  of X . The 2 energy is a wave is proportional to X  to  X   . Thus

Total energy = X 2 4π r  r2  Conservation of energy requires that this quantity remain constant as r  increases.  increases. Thus X2 α 

1 r 2

and

X α 

1 r 

As a seismic wave travels in the Earth the amplitude will decrease as 1/r , even if no energy is lost. This phenomena is known as spherical divergence or geometric spreading.

(b) Attenuation

As a wave passes through the Earth, the particle motion causes the material to be distorted and the wave energy is converted in heat. This results in an additional loss of energy, that is described  by an exponential decay:  X = X o e-kr 

16

Where e = 2.718, X  2.718,  X o is the amplitude at r=0 and r=0 and k is a constant. If k  is  is small, the attenuation will  be small, s mall, as k   increases, the attenuation becomes stronger. In a distance 1/k  1/k  the   the amplitude falls from Xo to Xo

1 e

.

Another common definition is the absorption coefficient, α, expressed in decibels per wavelength. This is based on the observation that the energy lost is dependent on the number of oscillations per second produced by the wave. Thus high frequencies will attenuation more quickly than low frequencies. This is illustrated in the MATLAB script waves_attenuation.m.

A consequence of frequency-dependent attenuation is that the shape of a pulse can change as it  propagates through the Earth. This occurs occurs for the following reasons: ●  The pulse can be decomposed into a set of sine or cosine curves with different wavelengths (Fourier’s Theorem and illustrated in MATLAB script fourier_square_wave.m) ● As the pulse travels the short wavelength signals attenuate more quickly. ●  The long wavelengths dominate, giving the pulse a smoother shape and longer duration (example in Kearey, Figure 3.7)

(c) Scattering

Suppose a medium is inhomogeneous and contains some grains with a different seismic velocity to the host rock. Seismic waves will be diffracted / scatteredfrom these grains and energy will be lost from the coherent wavefronts and turned into random seismic energy. The net result is that energy will be lost.  Footnote : Decibels A seismic wave changes in amplitude from A 1 to A2  as it travels from point 1 to point 2. The corresponding intensity changes from I 1 to I2. Note that I1 = A12. This change in decibels can be expressed as :

⎛  I   ⎞ ⎛  A  ⎞ dB = 10 log10 ⎜ 2 ⎟ = 20 log10 ⎜ 2 ⎟ ⎜ ⎟ ⎜ ⎟ ⎝  I 1  ⎠ ⎝  A1  ⎠

17

C1.8 Seismic energy sources

Basic requirements for an efficient seismic source include: ●Generate sufficient energy in an appropriate frequency band ● Economical  ● Non-destructive ● Repeatability

1.8.1 Land seismic sources

 Near surface exploration ● A high frequency seismic source is needed to generate short wavelength signals (the wavelength should be less than thickness of layers to be studied). However, if the frequency is too high, then attenuation will become a significant problem.

● A hammer striking a metal plate makes an  portable, repeatable source. A sledge hammer gives good energy, but relatively low frequency. A smaller hammer will give high frequency signal for shallow imaging. By stacking a number of hammer  blows, the signal-to-noise ratio can be increased.

●  Firearms give higher frequencies, but are not always permitted in civilized countries ( e.g. shotgun, semi-automatic etc.) ● Buffalo guns and Betsy guns allow a blank cartridge to be fired without using ammunition. ● Weight drops are also used and are commonly mounted on the back of a small truck. Larger weight drops have included the use of helicopters.

 Explosives Explosives were once widely used, but this has declined with the development of the Vibroseis method. Typically rotary drilling is used to place the shot in an 8-10 cm diameter hole at a depth of 6-30 m. Drills can be mounted on a truck, or slung beneath a large helicopter. While spectacular, blasting rocks and soil into the air represents a loss of energy and a safety hazard.  Advantages Sharp impulsive signal. Strong signals can be generated if needed for imaging the lower crust and upper mantle.

18

 Disadvantages  Disadvantages  Not repeatable, and quite slow. Most energy is in P-waves (why?). Many environmental concerns and cannot be used in urban areas .

 Drilling shot holes in Tibet, 1994 1994

1000 kg shot

Vibroseis

 Lithoprobe SNORCLE SNORCLE transect in B.C.

Truck generating signal 

For deeper exploration on land, the Vibroseis method is now the most widely used technique. This does not generate a sharp impulsive signal, but generates a longer waveform that typically lasts 10-20 seconds with frequencies 10-100 Hz. Data processing requires that the field recording is analysed for correlation with the known source signal. This can allow signal detection in areas with significant seismic noise (e.g. vehicle traffic, wind etc) Vibroseis trucks can generate forces up to 100,000 N and generally need a hard surface upon which to operate. It is debatable whether they cause no damage to a road. For extra signal strength, multiple trucks are used and they vibrate in phase. Both P-waves and S-waves can be generated.

19

Example 0 : Two separate reflections

Example 1 : Two overlapping reflections

20

Example 2 : Two overlapping reflections and moderate noise levels

Example 3 : Two overlapping overlapping reflections and and severe noise noise

21

Stacking is used to increase the signal-to-noise ratio. The inherent assumption of stacking is that the signal is the same each time ( coherent) and the noise is different ( incoherent). Coherent noise can not be removed by stacking.

22

1.8.2 Marine seismic sources

 Air guns An airgun works by releasing a bubble of high pressure air into the water. The rapid expansion of the bubble generates seismic energy with a frequency content around 10-100 Hz.

 Airgun array on Veritas ship

Operation of Bolt airgun. From Telford ch. 4. The chamber is filled with air at 10-15 MPa which is then released into the water.

 Airgun array and streamers streamers being towed behind a survey vessel.

 Advantages ● Very repeatable, reliable source.

 Disadvantages  Disadvantages

long wave train. However, by using an ●The bubble pulse oscillates, generating a relatively long array of air guns at differing depths, the combined waveform can be made shorter in duration (figure from Kearey). ● In water only P-waves are generated, but S-waves can be generated by mode conversion at the seafloor. ● In recent years there has been a lot of concern about how marines seismic exploration affects marine mammals.

23

 Explosives Once widely used in marine seismics, but now largely replaced by air guns. Bubble pulse problem is difficult to overcome (Kearey figure 3.15), but can be minimized by detonating near the surface (with some loss of energy). Occasionally large shots are used in deep studies but many restrictions are in place to protect oil production facilities and fisheries.  Higher frequency marine sources A range of higher frequency sources are used in marine exploration (100-10000 Hz). These include pingers, sparkers and boomers. The high frequency gives a sharp image of shallow sediment layers, but higher attenuation limits the depth of penetration. These methods are widely used in offshore geotechnical surveying (e.g. choosing locations for seafloor pipes and cables, in harbor construction).

1.8.3 Frequency content of seismic sources  Earthquake surface waves  Earthquake body waves Vibroseis  Air guns  Explosives

0.1-0.01 10-0.1 10-100 10-100 10-300

Hz Hz Hz Hz Hz

24

C1.9 Seismic detectors C 1.9.1 Onshore seismic exploration

Geophones On land, the surface moves as a P-wave or S-wave arrives. Generally reflected signals arrive at steep angles of incidence. Thus P-waves produce surface motion that is dominantly vertical. Geophones measure ground motion by converting motion into electrical signals. Most geophones measure a single component (vertical), but multiple component ones are sometimes used.

The most common design is the moving-coil geophone. Ensuring good mechanical coupling with the ground is essential, usually through a metal spike. Sometimes a cluster of geophones is planted and the output summed to improve signal to noise ratio. For shallow studies, geophone spacing can be as little as 1 metre. In deeper studies a spacing of 10-50 metres is more usual.

Geophones are manufactured to detect a particular frequency band. This should match the seismic source being used in a particular survey.

Geophones are connected to telemetry cables that transmit the recorded signals back to a recording unit.

Seismometers For earthquake studies a more permanent installation is usually required. Three components are usually recorded and the sensor is tuned to detect lower frequencies. Often the seismometer is  placed in a shallow vault to minimize wind and other forms forms of noise.

25

1.9.2 Offshore seismic exploration Hydrophones and streamers S-waves cannot travel in water, so only P-waves can be detected in the water. In marine exploration seismic waves are detected by the change in pressure as a P-wave passes the detector. This type of sensor is called a hydrophone  and it converts a change in water pressure into an electric signal through the piezoelectric effect.

In offshore exploration, many hydrophones are placed in neutrally buoyant, oil-filled streamers that can be in excess of 6 km long. The depth is controlled by movable fins and a tail buoy is used to determine the direction of the streamer relative to the ship.

Feathering occurs when ocean currents push the streamer at an angle to the direction in which the survey ship is travelling. Modern ships can pull multiple streamers, and this overcomes feathering and also provides 3-D coverage.

Ocean-bottom seismometers (OBS’s) These are like land seismometers, but are dropped to the seafloor from the sea surface. Coupling with the seafloor allows 3 components of motion to be recorded (i.e. P-waves and S-waves can be detected). Generally slow to use, and most widely used in research surveys. Some applications in hydrocarbon exploration.

Dalhousie University OBS on deck

Scripps Institution Institut ion of Oceanography OBS

Geophysics 224 MJU 2006

26

C2.1 Travel time curve for a single horizontal interface

Seismic energy can travel from the shot (*) to the receivers (geophones) by the 3 distinct routes shown in the upper panel. The lower panel is a travel-time curve and is a plot of travel time as a function of distance from the shot ( x).

(1) Direct wave A P-wave travels in a straight line just below the surface of the Earth. The travel time is t d  =  x / v1 which plots as a straight line on the travel time curve. This line passes through

the origin on the travel time curve ( x = 0 and t  =  = 0). (2) Ground roll A Rayleigh wave travels in a straight line along the surface at a velocity = v R. The travel time is t g =  x / v R which plots as a straight line on the travel time curve. This line passes

through the origin ( x = 0 and t  =  = 0). Since vR  < v 1 note that slower seismic signals have a steeper slope on this plot.

1

(3) Reflection

A P-wave reflects from the interface between layer 1 and layer 2. The angle of incidence and reflection are equal (see C1.5). Using Pythagoras’ Theorem, the distance travelled by the seismic signal on the downward leg of the journey is : d  =  z + 2

 x 2

4

From symmetry, the total distance travelled is 2d . The whole journey is travelled at velocity v1, so the travel time is given by  z + 2

t ref  = 2

 x

2

4

v1

=

4 z 2 + x 2 v1

Note that:

(1) When x = 0, t ref  ref  is not zero, in contrast to the direct wave and ground roll. (2) t ref  ref  has a minimum value when  x = 0. In this situation, the seismic signal travels vertically and makes an angle of 90° with the interface. This geometry is called normal incidence and the travel time is t ref  ref  = t 0 = 2 z / v1 (3) In practice, it is impossible to measure t 0  since a geophone will be destroyed if  place close to the shot. Observations at a range of  x  values can be used to extrapolate to find t ref  ref  at x = 0. The travel time can be written:  x 2 2 2 t ref  = t 0 + 2 v1 2

2 Thus a graph of t ref  versus  x  can be used for this. This will be a straight line with

intercept t 02 when x = 0 (4) At normal incidence, two unknown model parameters ( z and v1) determine t 0. The solution of this inverse problem (to find the model parameters  z and v1) is nonunique. This is because we have two unknowns and just one equation. (5) As x gets large the path taken by the reflection becomes very close to that taken by the direct wave. This means that the travel times for the reflection and direct wave  become very close and tref   x / v1 →

Can you sketch the travel-time curve for the reflection on the figure above? To address the problem of non-uniqueness in (3) and find the depth ( z) and velocity ( v1) we need to consider values of  x > 0.

2

2 ref 

t 

=

4 z 2 2

+

v1

 x

2 2

v1

= t  + 2 0

 x

2 2

v1

 Note that this equation is for a parabola (t ref  ref  varies as x squared). Simple re-arrangement gives : t ref  =

4 z 2 + x 2 v1

=

⎡ ⎣

2 z v1

t ref  = t 0 ⎢1 + (

1+ (

⎤ 2 z ⎦  x

)2 ⎥

 x

2 z

1/ 2

) 2 = t 0 1 + (

 x

2 z

)2

1

⎡  x 2 ⎤ 2 ) ⎥ = t 0 ⎢1 + ( v t  1 0 ⎣ ⎦

We can simplify the equation for t ref  ref   by using a power series expansion (Taylor’s Theorem) and assuming that x / v1t 0 is relatively small. 1

1

⎡ ⎡ ⎡ 1  x 2 ⎤  x 2 ⎤  x 2 ⎤ 2 = t 0 ⎢1 + ( ) ⎥ = t 0 ⎢1 + ( ) ⎥ = t 0 ⎢1 + ( ) + ...... ⎥ v1t 0 ⎦ v1t 0 ⎦ ⎣ ⎣ ⎣ 2 v1t 0 ⎦ 2

t ref 

If the higher order terms are ignored, then we can write that:

⎡ 1  x 2 ⎤  x 2 t ref  = t 0 ⎢1 + ( ) ⎥ = t 0 + 2 v t  2 2v1 t 0 1 0 ⎣ ⎦ Re-arranging gives an expression for t ref ref - t 0 which is termed the normal moveout. t ref  − t 0 =

 x 2

2v12 t 0

 Normal moveout (NMO) is a measure of the extra time taken for seismic signal to travel on a non-vertical path, compared to the time for a signal travelling vertically. The word normal refers to the seismic energy travelling at right angles to the interface. 2

A graph of NMO versus  x  will be a straight line, passing through the origin. The slope of the line will be 1 / 2v12 t 0 . As  x increases, the approximation made above becomes invalid

and a deviation from the straight line will be observed.

Example 1

In this example, we will look at the synthetic seismic reflection data and compute the depth of the interface and estimate the velocity of the layer. These data are shown in graphical and numerical form on pages 4 and 5. (1) Read the value of t 0 from the table

t 0 =

________________

(2) Look at the geophone at  x = 90 m.

tref  =

________________ ________________

3

(3) At this geophone, compute the NMO

NMO = ________________ ___ _____________

(4) Find v1 using the equation above

v1 =

________________ ________________

(5) Using the values for v1 and t 0 find the depth of the interface. This can be done with the equation for travel time at normal incidence t 0 = 2 z / v1 . Rearrange this equation to find the depth to the interface.  z = __________________ Does your answer agree with that shown in the plot below? If not, can you suggest a reason for the discrepancy? (6) What can we determine from the seismic data about the seismic velocity below the interface (v2)?

4

2

Note : To be more rigorous when computing v1 it is best to plot a graph of  NMO vs.  x which should be a straight line with slope 1 / 2v12t 0 . A non-straight line may indicate that

the simple form of the NMO equation is invalid. C2. 1 1- r ef l x( m) - 39 390. 0. 00 000 0 - 37 375. 5. 00 000 0 - 36 360. 0. 00 000 0 - 34 345. 5. 00 000 0 - 33 330. 0. 00 000 0 - 31 315. 5. 00 000 0 - 30 300. 0. 00 000 0 - 28 285. 5. 00 000 0 - 27 270. 0. 00 000 0 - 25 255. 5. 00 000 0 - 24 240. 0. 00 000 0 - 22 225. 5. 00 000 0 - 21 210. 0. 00 000 0 - 19 195. 5. 00 000 0 - 18 180. 0. 00 000 0 - 16 165. 5. 00 000 0 - 15 150. 0. 00 000 0 - 13 135. 5. 00 000 0 - 12 120. 0. 00 000 0 - 10 105. 5. 00 000 0 - 90. 00 000 0 - 75. 00 000 0 - 60. 00 000 0 - 45. 00 000 0 - 30. 00 000 0 - 15. 00 000 0 0. 00 000 0 15. 15. 000 000 30. 30. 000 000 45. 45. 000 000 60. 60. 000 000 75. 75. 000 000 90. 90. 000 000 105. 105. 000 000 120. 120. 000 000 135. 135. 000 000 150. 150. 000 000 165. 165. 000 000 180. 180. 000 000 195. 195. 000 000 210. 210. 000 000 225. 225. 000 000 240. 240. 000 000 255. 255. 000 000 270. 270. 000 000 285. 285. 000 000 300. 300. 000 000 315. 315. 000 000 360. 360. 000 000 375. 375. 000 000 390. 390. 000 000

ect i on, exam xampl e 1 r ef 1( ms) 593. 59 3. 33 337 7 589. 58 9. 02 021 1 584. 58 4. 85 856 6 580. 58 0. 82 822 2 576. 57 6. 92 922 2 573. 57 3. 20 200 0 569. 56 9. 59 596 6 566. 56 6. 17 172 2 562. 56 2. 89 891 1 559. 55 9. 77 773 3 556. 55 6. 81 819 9 554. 55 4. 02 029 9 551. 55 1. 40 407 7 548. 54 8. 95 951 1 546. 54 6. 66 666 6 544. 54 4. 56 561 1 542. 54 2. 62 628 8 540. 54 0. 87 876 6 539. 53 9. 29 299 9 537. 53 7. 91 911 1 536. 536. 69 699 9 535. 535. 67 671 1 534. 534. 83 832 2 534. 534. 17 177 7 533. 533. 70 709 9 533. 533. 42 427 7 533. 533. 33 333 3 533. 533. 427 427 533. 533. 707 707 534. 534. 177 177 534. 534. 832 832 535. 535. 671 671 536. 536. 699 699 537 537.. 905 905 539 539.. 299 299 540 540.. 868 868 542 542.. 628 628 544 544.. 561 561 546 546.. 666 666 548 548.. 951 951 551 551.. 394 394 554 554.. 029 029 556 556.. 819 819 559 559.. 773 773 562 562.. 891 891 566 566.. 172 172 569 569.. 596 596 573 573.. 180 180 584 584.. 833 833 589 589.. 021 021 593 593.. 337 337

5

(1) Look at the seismic reflection data presented above. How are they different from those in Example 1? (2) Which model parameters ( v1, v2 and z) will be different?

6

7

C2.2 Travel time curves for multiple reflectors ●

The ray paths and travel times for multiple layers can be computed using ray-tracing, as demonstrated in Lab 3.

●

MATLAB script reflect_layers_v1.m performs simple ray tracing.

x(m) ref1(ms) ref2(ms) ref3(ms) 0.000 533.333 933.333 1219.048 30.000 533.707 533.707 933.455 1219.110 60.000 534.829 534.829 933.819 1219.298 90.000 536.695 536.695 934.426 1219.610 120.000 539.296 935.276 1220.047 150.000 542.623 936.367 1220.608 180.000 546.659 937.698 1221.294 210.000 551.394 939.269 1222.103 240.000 556.807 941.074 1223.037 270.000 562.881 943.119 1224.094 300.000 569.589 945.395 1225.274 330.000 576.918 947.904 1226.575 360.000 584.833 950.642 1228.000 390.000 593.322 953.609 1229.545

1

●

The two examples both both show three reflections. Note that: (a) Rays are refracted at interfaces according to Snell’s Law (b) The polarity and amplitude of a reflection is due to the impedance change at the interface (c) Deeper reflections are flatter (i.e. (i.e. they  they exhibit less normal moveout) (d) The high amplitude reflection from a shallow interface can weaken later reflections (from deeper structure) (e)For a shot gather, the ray paths reflect (bounce) at different points on the reflector surface. (f)This geometry, with geophones on each side of the shot is called a split-spread

2

●

In the real shot gather gather note that:

(a) The vertical scale is two-way travel time in seconds. (b) Positive arrivals have the trace shaded black. This helps the eye to correlate traces when noise is present. Beer sometime helps with this process, but too much has  been known to have the reverse effect. Calibration is needed to optimize this task. (c)  Note the direct wave, which is the first arrival. (d) The ground roll is long lasting and obscures reflections at small offset (A). (e) Slower arrivals plot at a steeper angle on the seismogram (more time needed to travel the same distance) (f) Reflections can be recognized by their parabolic shape and are labelled B,C,D, and E. Note that they are never the first arrival. This means that reflections are recorded when the ground is still moving from earlier arrivals, and reflections can  be obscured by noise.

3

Interval velocity and average velocity

Consider the multiplayer Earth model shown above. The ith layer has a thickness  z i and a velocity vi . The ray travels through the layers 1 to n and is refracted at each interface. It is then reflected by the nth interface. The seismic wave spends a time τi in the ith layer. The interval velocity is the actual velocity in a specific layer and is defined as vi = z i / τ i For the whole ray path we can define an  average velocity as

V  =

∑ z  = ∑ v τ  ∑τ  ∑τ  i

i

i

i

i

Another way of averaging the velocity is to use the root-mean-square average. This is defined as: 1

⎡ ⎤2 2 τ  v ⎢∑ i i ⎥ ⎥ = ⎢ i =1n ⎢ ⎥ τ i ∑ ⎢⎣ i =1 ⎥⎦ n

V rms ,n

and is needed to compute interval velocities, as described below. For the case of n = 2

V rms, 2

⎡v =⎢ ⎣

1

+v τ 1 + τ 2

2 1 τ 1

2 2 τ 2

⎤2 ⎥ ⎦

4

Normal move-out for multiple layers and the Dix equation

To determine the depth of the interfaces and the interval velocities, an approach similar to that used in C2.1 is needed. Note that we need more than just the travel-time at normal incidence, and should use a plot of NMO vs. offset ( x).  x). When the seismic signals travel close to the vertical direction, can show that the normal moveout for the nth reflection is:

Δt  = t − t n =

 x

2

2 2t nV rms ,n

This is very similar to the equation in C2.1 for the normal moveout above a single interface. The only change is that the interval velocity (v 1) has been replaced by the r.m.s. velocity (i.e. the average velocity down to the nth reflector). It can be shown the rootmean-square velocity, gives the correct average in this context.

Consider the seismogram with three reflections that is shown above: For each reflection we can measure the zero-offset travel times (t   (t 1  , t 2 , t 3). If necessary a 2 2  plot of t   versus x  versus x  may be needed (see C2.1)  Next the r.m.s. velocities  are computed from the normal moveout. Consider the 4 th geophone at distance x distance  x.. The normal moveout at this value of  x are  x are  Δt 1,  Δt 2 and  Δt 3. From these NMO values, the average (r.m.s.) velocities can be computed.  x 2  NMO = Δt 1 ( x) = t ( x ) − t 1 = 2 2t 1V rms ,1 Re-arranging this equation gives V rms,1 =

 x 2 2t 1 Δt 1 ( x)

5

For the other reflections it is obvious that  x 2

V rms, 2 =

2t 2 Δt 2 ( x )

and V rms ,3 =

 x 2 2t 3 Δt 3 ( x )

The next stage of the data analysis is to compute the interval (true) velocities (v ( v1 , v2 , v3) from the average (r.m.s.) velocities (V  (V rms,1 rms,1 , V  rms,2  rms,2 , V  rms,3  rms,3) The is achieved with the Dix equation that states:

⎡V  vn = ⎢ ⎢⎣

1

t  ⎤ 2

t  − V 

2 rms , n n

2 rms , n −1 n −1

t n − t n−1

⎥ ⎥⎦

While this appears complicated, it is really a simple recursion equation. This means that if we know the velocity in one layer, the equation tells us what the velocity will be in the layer below. To start this process, we begin at the surface. For the first reflection, the seismic signal travels only in Layer 1. Thus v1= V rms,1 rms,1  Now V rms,2 rms,2 represents an average velocity of layers 1 and 2. Thus to find v2 we can use the Dix equation with n = 2

⎡V  v2 = ⎢ ⎢⎣

1

t  − V 

2 rms , 2 2

t  ⎤ 2

2 rms ,1 1

t 2 − t 1

⎥ ⎥⎦

This can then be repeated as many times as necessary. For the next layer down we can write

⎡V  v3 = ⎢ ⎢⎣

1

t  − V 

2 rms, 3 3

t  ⎤ 2

2 rms , 2 2

t 3 − t 2

⎥ ⎥⎦

Example  : Look at the data in the table above for the three reflection seismogram in Example 1. The travel times for the first reflection are the same as for the example in C2.1

The solution to that problem gave:

v1 = 1500 m s -1 ;  z 1 = 400 m and t 0 = 0.533 s -1  Note that V rms,1 rms,1 = v1 = 1500 m s

6

Consider the second reflection Choose a geophone,

 x

= _____________ _____________

Zero offset travel time,

t 2

= _____________ _____________

 Δt 2 

=______________

V rms,2 rms,2

=______________ 

 Normal moveout r.m.s. velocity for Layers 1 and 2

Use the Dix equation to compute the interval velocity for layer 2: 1

⎡V  v2 = ⎢ ⎣⎢

t  − V 

2 rms , 2 2

t  ⎤ 2

2 rms ,1 1

t 2 − t 1

⎥   ⎦⎥

=_______________

Does this agree with the true model?  Now compute the thickness of the second layer. To do this note that the reflection spends a time t 2-t 1 in the second layer. In this time, it travels a distance 2z2 . Thus we can write

v2 =

2 z 2 t 2 − t 1

Re-arranging the equation gives

 z 2 =

v2 (t 2 − t 1 ) 2

= _________________ _________________

 Repeat this for the third reflection reflection Choose a geophone,

 x

= ______________ ______________

Zero offset travel time,

t 3

= ______________ ______________

 Δt 3 

=______________

V rms,3 rms,3

=______________ 

 Normal moveout

⎡V 

t  − V 

⎢⎣

t 3 − t 2

v3 = ⎢

 z 3 =

1

2 rms, 3 3

v3 (t 3 − t 2 ) 2

t  ⎤ 2

2 rms , 2 2

⎥ ⎥⎦

= _______________ _______________

= _________________ _________________

7

Multiple reflections

(a) The second reflection is a multiple bounce from the interface at 0.4 km depth. (b) The travel time of the multiple at zero-offset is exactly double that of the single reflection. This can help identify a multiple. (c)  Note that when the seismic energy bounces at the Earth-Air interface, the reflection coefficient can be large. Remember that

 R =

−  ρ 1v1  ρ 2 v 2 +  ρ 1v1  ρ 2 v 2

Consider the upward travelling wave as it approaches the surface. ρ2 is now the density of the air which is essentially zero. This gives R = -1 and almost all energy is reflected downwards. (d) For an example of multiples in single channel seismic data, see Kearey Figure 4.56 and the discussion in C2.3

8

C2.3 Strategies for collecting seismic data C2.3.1 Shot gather

The shot gather allows us to determine the variation of velocity with depth through the measurement of the normal moveout.

●

●

The seismic survey measures time. Processing techniques are used to convert time to depth.

 In a shot gather, each ray will reflect at a different location. If the structure varies significantly in the horizontal direction, then this will introduce errors.

●

C2.3.2 Single channel seismic profiling

•

•

•

In many seismic surveys, we are most interested in how structure varies with horizontal position and may be less interested in actual depths.  Note that the seismic section is composed of many individual traces . A trace is the time variation of the ground motion (or water pressure in a marine survey) recorded by the receiver for a single shot. In this type of seismic survey the absolute depth cannot be determined. The travel time = t = depth/ velocity. With just one measurement of travel time, t, it is

1

impossible to determine two unknowns. Only relative changes in depth can be determined. •

•

Seismic reflections are generally quite weak (low amplitude) and this technique does not generally work well on land where cultural noise levels are higher and the coupling between the geophone and the ground is variable.

Example above is from marine seismic exploration in the Gulf of Corinth (Kearey Figure 4.56). Note that D is a direct wave where a signal travelling horizontally  passes the receiver. All other arrivals have travelled vertically downwards in the Earth and back again. The seabed (SB) gives a strong reflection, and multiple reflections are quite strong (SBM1 and SBM2) and obscure other reflections from greater depth. RH is a reflection from bedrock that is exposed above the seafloor on the right side of thee profile.

2

C2.3.3 Multi-channel seismic profiling

•

•

•

By combining the two techniques described above, we can measure both the velocity variation with depth and   horizontal variations in structure   on a  profile. One simple acquisition strategy is sketched above. Note that at location ‘A’ the reflector is quite shallow. At location ‘B’ the reflector is deeper and the reflector shows a larger zero-offset travel time and less normal moveout. Analysis of these data will give estimates of both the interface depth and velocity. However this has the disadvantage that the energy received at each receiver has  been reflected at different points. The common mid-point technique can overcome this problem. Seismic data is grouped into shot-receiver paths that share a common reflection point. This technique is discussed in detail in C2.5

 MJU Geophysics 224 2006 

3

C2.4 Travel time curve for a dipping reflector

t  =

 x 2

+

4 z 2

+

4 xz sin θ 

v1

●

t(x) is greater than t(-x) and the travel time curve is asymmetric about x = 0

●

the minimum travel time does not occur at x = 0 m (why?)

●

also note that the reflection received at x = 0 m did not originate beneath  x = 0.

●

To account for this effect a technique called migration is used. More on this later! Using the same approach as in C2.1, we can show that

t  ≈ t 0

+

( x 2

+

4 xz sin θ )

2v12 t 0

The dip moveout is defined as

ΔT d  =

t  x

− t − x =

2 x sin θ  v1

 MJU Geophysics 224 2006 

1

C2.5 Common mid-point profiling ● The

shot gather has the disadvantage that each ray has been reflected from a different  point on the interface. If there are lateral variations in structure, this will result in errors.

● This

can be overcome by using multiple shot points and multiple receivers. From the total dataset, a set of rays are then chosen that have a common reflection point.

● This

technique is called both common-mid point (CMP) or common depth point (CDP)  profiling. The shot and receiver are located equal distances from the CMP.

Simple example of a CMP gather

Matlab script cmp_v1.m shows the path geometry for a CMP gather.  Note that the offset is the distance from shot (*) to receiver (v). Thus with geophones from –200 m to + 200 m, the maximum offset between shot and receiver is 400 m.

●

The parabolic travel time curve derived for the shot gather can also be used to analyse a CMP gather.

●

1

The technique for collecting CMP data is illustrated above. A seismic source is located located at the front of a 6-channel array.

●

The spacing of geophones is

●

Δx

along the entire array.

After each shot is fired, the array is moved forward a distance n Δx.

●

In a marine survey this would consist of the whole airgun and hydrophone streamer  being towed forward by the survey ship. On land the last geophone would be brought up to the front of the array, a new hole drilled, loaded with explosives and fired.

●

● If

we have N geophones and that array moves a distance n  Δ x between shots, shots, then you can show that the number of rays that share the same common mid-points  = N/2n. This N/2n. This quantity is also called the fold or the coverage (in percent).

● This

survey will give 3 rays for each mid-point. This is called 3-fold CMP coverage  or 300% coverage.

● Kearey

Figure 4.14 show an example with 6-fold CMP coverage .

2

Stacking ●  Once

we have collected a number of rays that share the same mid-point, we can use stacking to improve signal to noise ratio in the data. This is needed because reflections are generally weak. Remember that stacking adds a number of traces. If the signal is the same in each, but the (random) noise is different, then the noise will cancel and the signal be enhanced.

● Prior

to stacking, the individual traces must be shifted vertically (i.e. in time) to allow the reflection to occur at the same time.

● Traces

at non-zero offset are shifted in time to be equivalent to the zero offset trace.

●  This

requires that we compute the expected NMO for a given reflection and offset. However we don’t know the velocity (this is the goal of the survey!) so a trial and error approach is used.

●  This

is illustrated in cmp_v2.m. If the stacking velocity is too high, then the curve remains a downward parabola. If the stacking velocity is too low, then the curve is an upwards parabola.

● Each

reflection is analysed and the stacking velocity systematically varied. The velocity that gives the best stack (flattest reflection) is then taken to be the r.m.s. velocity from the surface to the reflector. This is equivalent equivalent to plotting NMO NMO as a function of of x2.

● A

different velocity will be needed for each reflector.

● Note

that the direct wave and ground roll do not stack coherently, and this process tends to remove them from the data.

3

4

5

6

Geophysics 224 MJU 2006 

7

C2.6 Resolution in seismic reflection surveys

When interpreting geophysical data, it is very important to be aware of the resolution that is possible with the method being used.

C2.6.1 Vertical resolution

Consider the thin layer shown in the Figure C2.6 with the upper surface at a depth of 0.2 km. Seismic reflection data is collected at normal incidence. The reflections from the upper and lower surfaces were computed with MATLAB script vert_resn_reflection.m

t  =  = 200 m The seismic pulse has a wavelength, λ =200 =200 m. When the layer is relatively thick, the two reflections are distinct in the combined seismic trace (right hand  panel).

t  =  = 20 m When the layer is very thin (t  (t   = 20 m), the two reflections overlap and two separate reflections cannot be identified. In this case, it is not  possible to determine that two layers are  present.

1

t  =  = 50 m  Note that the transition occurs when the crest of second arrival coincides with the trough of the first arrival. This is produced by the second reflection travelling a distance 2t  greater than the first reflection.

t = 60 m The second reflection can just be detected.

Two waves can be distinguished from each other if they are separated by distance  λ /2.  Now the reflection from the lower interface has travelled an additional distance of 2t . Thus the best resolution occurs when 2t  =  =  λ /2 which  /2 which gives t =  λ /4

Example: A 40 Hz seismic signal reflects from the top and bottom of a layer in a medium with v ~ 3000 m/s. What is the minimum layer thickness that could be detected?

2

How does this compare with features seen in an outcrop?

C2.6.2 Horizontal resolution

Two factors can limit horizontal resolution in a reflection survey. One is due to the  physics of seismic wave propagation (Fresnel zones), and the other is due to the geophone spacing. The resolution will be which ever number is larger.

Fresnel zones

It is sometime tempting to imagine that a seismic reflection occurs only at the point where the ray hits an interface. However it is important to remember that a complte understanding of seismic wave propagation requires us to consider both wavefronts and rays. Consider seismic energy that strikes the interface at a distance from the  point below the source (where normal reflection occurs). These signals are not predicted to return to the source by ray theory, but are expected according to Huyghens Principle. These arrivals will contribute to the seismic energy received at the source. However, they will travel a greater distance  than the normal reflection. Following the logic of the previous section, these arrivals will be in phase with the normal reflection until their round trip path length is λ /2 /2 greater than the normal reflection.

It can shown that the reflection essentially occurs from a circular disk, over which in  phase reflections occur. This is called the Fresnel zone  and has a width w

=

2 z λ 

3

The width of this zone is the smallest feature that can be detected with the surface. Note again that higher resolution will be achieved by using short wavelengths. However this will result in more attenuation, and the signals will have limited penetration. Example: A 40 Hz seismic signal reflects from the base of a 2000 m thick layer with v ~ 3000 m/s. What is the minimum size of horizontal feature that could be detected with this configuration?

Detector spacing

The reflection points on an interface will have a horizontal spacing equal to half the geophone/hydrophone spacing. By making this spacing less than  the Fresnel zone, the survey resolution will not be limited by the layout, but by the physics of the wave  propagation.

Geophysics 224 MJU 2006 

4

C2.7 Statics, noise reduction and filtering of seismic reflection data C2.7.1 Statics and their removal

So far we have assumed that onshore seismic data is being collected on a flat surface with uniform near surface geology. However this assumption is not often valid, and two effects will cause a variation in arrival times between closely spaced traces. These time delays are called statics because they delay (or advance) all events in the seismic trace recorded by a particular shot-geophone combination.  Elevation statics:  geophones at an elevation above a datum (or reference level) will detect the incoming signal later than a geophone on the datum.

Weathering statics : statics : weathering of near surface rocks produces a zone of low velocities that is variable in thickness

Kearey Figure 4.15

1

Why are statics a problem in seismic reflection data?

Statics produce a vertical offset in seismic traces, as shown in the CMP section generated by cmp_v3.m. Thus when the traces are stacked, they will not stack coherently.

2

3

 Example 1 Statics are present in 50 Hz seismic reflection data. How large can the statics be and still permit coherent CMP stacking? Signal period (T) = ___ ms For reinforcement, maximum time delay for reinforcement = T/4 = ____ ms

 Example 2 If there is no weathered layer, and the bedrock has v = 3000 m/s, what elevation change will produce a static of 5 ms?

 Example 3 The weathered layer has v=1000 m/s and the bedrock has v = 2500 m/s. What change in the thickness of the weathered layer will produce a delay of 5 ms?

Correction of statics ●Elevation

statics can be removed by accurate surveying with GPS or differential GPS. A vertical accuracy of ±1 m is usually sufficient.

●Weathering

statics can be estimated from a variety of techniques, including measuring the refracted arrival (head wave) that travels along the base of the weathered layer. Details in Kearey, chapter 4, 58-59 .

● In

residual static analysis, traces are automatically aligned to produce the most continuous seismic event. Example below is from Kearey Figure 4.16

Before

After residual static analysis

4

C2.7.2 Frequency domain filtering

When seismic reflection data is recorded, the individual traces data are contaminated by signals other than seismic reflections. These signals include ground motion due to wind noise, vehicles etc. There are also other seismic phases that can obscure the reflections (ground roll). Electrical powerline noise (50 Hz and 60 Hz and harmonics) can also be  picked up by the telemetry cables. In marine surveys, noise signals originate in the motion of the ship and streamer through the water. Some of this noise has a frequency content that is quite distinct from the seismic signals  being recorded. These can be removed by frequency domain filtering techniques. Notch filter : removes a very narrow band of signals (e.g. 60 Hz powerline noise) High pass filter: All frequencies above the filter frequency are transmitted. This an remove ground roll and certain ship generated noise (also called a low-cut filter) Low pass filter: all frequencies below the filter frequency are transmitted. Can be used to remove high frequency noise in some circumstances. Band pass filter : transmits signals in a narrow frequency band.

 Example:  Example: Kearey figure 4.19. Note that the frequency content in the seismic section changes with time. Low frequencies dominate at later times, as the high frequencies are more rapidly attenuated than lower frequencies.

However, many artefacts in a seismic section have a frequency content that is very close to those present in the reflections. Thus they cannot be removed by conventional frequency domain filtering.

5

C2.7.3 Convolution

Convolution is a mathematical operation that can be defined as  y(t)

= f(t)*x(t)

=

∫  f (

) x(τ   − t )d τ  

τ  

Graphically, convolution is illustrated in the attached figure, where the function  x(t) is a series of spikes (delta functions) and  f(t)  is a sinusoidal wavelet (MATLAB script convolution_C27.m).

This mathematical procedure may seem quite abstract, but it occurs in reflection seismology, as illustrated in Kearey Figure 4.6 .

6

The sharp changes in acoustic impedance (Z) in the geological section are called the reflectivity function. The locations where this function has spikes are the locations that reflect seismic energy back to the surface. However, to measure the reflectivity function we need to send a seismic signal through the Earth. This seismic pulse has a finite length (in both space and time) and the seismic trace (measured by the geophone) can be shown to be the convolution of the reflectivity function and the seismic pulse . This has the effect to spread out the sharp response in the reflectivity function associated with each interface. As a consequence, the measured seismic trace is lower in resolution and we cannot accurately locate the time at which a particular reflection arrives. This will limit resolution, since we showed in section C2.6 that layers can only be resolved if they are separated by more than λ /4. /4. Thus a sharper response will be obtained by making the seismic pulse as short as possible in time. However, this requires an increase in frequency, and this will have the side effect of increasing attenuation. The mathematical process of deconvolution  is used to reverse the process described above. The goal is to find the sharp reflectivity function from the measured seismic trace.

7

C2.7.4 Deconvolution

The process of deconvolution can be used to undo the effect of other filtering operations. This includes the effect of a finite length seismic pulse, as illustrated above. Often used to shorten the length of the seismic pulse. A seismic pulse can be lengthened  by a number of processes: Reverberation : multiple reflections in the water layer. Ghosting : reverberation in the (low velocity) weathered layer at the surface.

To overcome these, and other effects, a range of deconvolution techniques may be used. More details and examples can be found in Kearey (p. 62-65). Kearey 4.22:  Reverberation removal with predictive deconvolution.

Before ●

After

remember that in Vibroseis ® surveys, the seismic pulse is very long. However, the fact that the shape of the pulse is known allows us to correlate the seismic trace with the known waveform during processing. However, in surveys using explosives are airguns, the source function is generally not well enough known for this to be applied.

Geophysics 224 MJU 2006 

8

C2.8 Migration of seismic reflection data C2.8.1 Why is migration needed? C2.8.1.1 Dipping layers are imaged imaged with an incorrect incorrect dip

Consider a shot (S) and receiver (R) that are located close together. A reflector has a true dip angle = αt  and is at a depth depth z   z  below   below the shot point. Seismic energy that returns to R will reflect at ‘A’ where the ray path is at 90º to the reflector. However, in a seismic section, a reflection in plotted as it was directly below  the shotreceiver point , which which in this case is ‘B’. Note the line A-B is the arc of a circle, centered at RS. The result is the reflector is imaged by the seismic data with an apparent dip of α s that is less than the true dip. From triangle A-C-RS  A-C-RS 

cos αt  = p / z

From triangle O-RS-C  O-RS-C 

tan αt  = z / x

From triangle O-RS-B  O-RS-B 

tan α s = p / x

Rearranging these equations gives

 sin αt  = tan α s

True dip (αt) 5º 10 º 20 º 40 º 60 º 80 º

Apparent dip ( α s) 4.98 º 9.85 º 18.88 º 32.73 º 40.89 º 44.56 º

The result is that reflectors appear to have an apparent dip that is less than the true dip . The effect is small at shallow angles, but becomes significant when the true dips are greater than 20 º. Note that a vertical reflector will appear at 45 º on a seismic section.

1

C2.8.1.2 Multiple reflections in same zero offset trace ●

  When a syncline has very limited topography, only one reflection is observed in each zero offset trace (i.e. shot and receiver are placed very close together and moved along the profile). As shown in the previous section, the dip will not be measured correctly. However, the overall geometry in the seismic record section will be very similar to the reflector geometry. Note that reflections occur at right angles to the surface.

  With more rugged topography, a normal reflection will occur at three locations when the source and receiver are above the syncline. This produces the characteristic “bow “bow tie” tie” in the travel times.

●

2

3

C2.8.1.3 Diffractions

•

Diffractions occur from point reflectors and corners. The figure above shows a small sphere that is diffracting seismic waves. The diffractor has the property that it scatters energy in all directions. Some of this energy returns to the location of the shot (*) along the same path (like a boomerang).

•

When the shot-receiver array is moved along the profile with zero-offset, the travel time curve above has a travel time t diff 

=

2  x 2

+

v1

h2

  which plots as a

hyperbola. •

A hyperbola will also result when the shot-receiver offset is varied (e.g. (e.g.   a shot gather or CMP gather). In this case the move out rate will be half that observed from a plane reflector

4

C2.8.2 Migration methods

Migration in geophysics has nothing to do with birds flying south in winter. It is mathematical process that attempts to reconstruct the reflector geometry from measurements of the reflected seismic energy. A wide range of methods have been developed.

The simplest technique works by assuming that from the travel time, we know that the reflection point must lie somewhere on a hemispherical surface. By migrating several  points this allows us reconstruct the reflector surface. If applied to a diffraction, this will collapse the hyperbola to a point . The diffraction hyperbola is the steepest dipping even that can be recorded on a seismic record section. Some commonly used methods of migration include: incl ude:  diffraction migration

●

●

 wave equation migration

 Kirchoff migration

●

More details and references for these methods are given by Kearey in chapter 4. Each migration technique can be further divided into time and depth migration time migration: migration: the vertical axis on the migrated section is still time.

●

depth migration: migration: estimates of velocity are used to convert time to depth during migration. Even rough estimates of velocity (e.g. ( e.g.   from computing the stacking velocity for each CDP gather) can greatly improve the quality of the migration.

●

Migration can also be applied before or after stacking the CMP gathers:  Post-stack depth migration  migration  : This was the original migration strategy. The zero offset traces derived from the CMP stacks are migrated. This has the advantage of lower computer time  and the higher signal-to-noise ratio can make the migration more stable. However the process of stacking assumes a relatively flat interfaces. This is valid in geological structures that have sub-horizontal stratigraphy, but is not valid in complex, highly 3-D environments.

●

●

 Pre-stack depth migration: This migrates individual traces prior to stacking. Requires major amount of computation, but is necessary in complex 3-D environments.

5

C2.8.3 Examples of migration

Wave equation migration in a location where structures generally subhorizontal. Note that bow-ties are essentially removed by the migration (Kearey chapter 4).

The Sigsbee model was developed to simulate sub-salt seismic imaging in deeper water exploration areas such as the Gulf of Mexico. More details can be found in Integrated approach to sub-salt depth imaging: synthetic case study, The Leading Edge, Edge , 21, 12171223, December 2002.

In the zero-offset time section, note the diffraction hyperbolae   that originate in the sharp corners in the upper surface of the salt sheet. Very few coherent reflectors can be seen beneath the salt.

6

Left : Kirchoff depth migration. Reflectors to the left and above the salt body are well imaged.  Right  : Finite difference Pre-stack depth migration. Coherent reflectors are imaged in the sedimentary sequence beneath the salt body.

Multiples do not migrate correctly and along with other errors in migration this can  produce “ smiles”  smiles” which are upward directed hyperbolae. This is illustrated in Figure 4.84 from Telford.

Original time section

Results of time migration

Example below is from a PGS advert in The Leading Edge. Note that structures with steeper dip are imaged more reliably with pre-stack depth migration.

Post-stack depth migration

Pre-stack depth migration Geophysics 224 MJU 2006 

7

C2.9 Hydrocarbon exploration with seismic reflection

Seismic processing includes the steps discussed in class: (1) Filtering of raw data (2) Selecting traces for CMP gathers (3) Static corrections (4) Velocity analysis (5) NMO/DMO (5)  NMO/DMO corrections (6) CMP stacking. (7) Deconvolution and filtering of stacked zero offset traces (8) Migration (depth or time)

Question : This sequence is for post-stack  depth migration. How will it change for prestack  depth  depth migration?

Processed seismic data can contribute to hydrocarbon exploration in several ways: ●Seismic

data can give direct evidence of the presence of hydrocarbons ( e.g.  bright spots, oil-water contact, amplitude-versus-offset anomalies).

●Potential

hydrocarbon traps  can be imaged (e.g.  reefs, unconformities, structural traps, stratiagraphic traps etc)

●Regional

structure can be understood in terms of depositional history and the timing of regression and transgression (seismic stratigraphy, seismic facies analysis). This can sometimes give an understanding of where potential source rocks are located and the relative age of reservoir rocks.

1

C2.9.1 Direct indicators of hydrocarbons 2.9.1.1 Bright spots As shown in C1.3, a gas reservoir can have a P-wave velocity that is significantly lower that the surrounding rocks. This can lead to a high amplitude (negative polarity) reflection from the top of the reservoir.

● However,

not all bright spots are hydrocarbons. They can be caused by sills of igneous rocks or other lithological contrasts. In areas of active tectonics, ponded partial melt can  produce a bright spot (e.g. Southern Tibet).

● It

should also be noted that true amplitudes are not always preserved in seismic data recording and processing. When data is recorded, the amplification of the signal varies from trace to trace. It can also vary with time in a given trace (later signals have travelled deeper and have a smaller amplitude so the amplification is increased). Thus variations from trace to trace in a seismic section do not necessarily imply a change in a sub-surface property.

●Focussing

and defocusing   effects on an undulating interface can also change the measured amplitude. Beware!

2.9.1.2 Hydrocarbon-water interface i nterface Example 1 : Fulmar field, North Sea

The oil-water contact is flat in this depth section from the Fulmar field in the  North Sea (Kearey 4-39). Note that it crosses reflectors in the anticline.

Example 2 : The Alba Field, North Sea

Details in Macleod et al., (1999)

2

Example 3 : Alexis Field, Alberta

●The

Alexis Field, 65 km NW of Edmonton, Alberta is described in Chapter 6 of  Andersen et al, 1989. The oil and gas reservoir is located in a structural remnant of dolomitized carbonate of the Mississippian age Banff formation. Upper surface is erosional in nature and overlying rocks are lower velocity Cretaceous shales and sandstones.

●The

seismic data use the reflection from the Viking Fm as a datum. This requires all the seismic traces to be aligned on this reflector. This removes any residual statics and the velocity effects of shallower structure. This is necessary because variations in the reflection from the reservoir are quite subtle.

●To

interpret the seismic reflections, velocity measurements are made in wells (sonic logs). From the variation of velocity with depth, a synthetic seismogram   can be computed. This is the predicted seismic trace that would be measured at this location. It allows us to determine which velocity and lithology changes will be detectable in surface seismic data.

●Remember

that the finite wavelengths of seismic data do not allow features less than λ /4 /4 in thickness to be imaged (section C2.6).

●The

reservoir in this case is quite thin, and completely separate reflectors are not seen from the upper and lower interfaces. However, the shape of the reflector changes distinctly between a wells in the reservoir (15-36) and a well outside the reservoir (3-1).

3

2.9.1.3 Amplitude versus offset (AVO)

Previously we computed reflection coefficients for seismic waves incident at normal incidence. However, the full Zoeppritz equations predict that the reflection coefficient will vary with the angle of incidence . The two examples below were generated with MATLAB script AVO_simple.m and illustrate this phenomena. Example 1

Both P-wave and S-wave velocities increase across the interface. Remember that at nonnormal incidence, an incident P-wave will generate 4 new waves. These are the reflected and transmitted P-waves and  the reflected and transmitted S-waves. The reflection coefficient of the P-wave decreases with angle in this example.

4

Example 2

This example simulates the effect of shale overlying a gas saturated sand reservoir. The Pwave velocity decreases , giving a negative reflection coefficient at normal incidence.

The S-wave velocity increases from the the shale into the gas sand. Why? Can show that this corresponds to a decrease  of Poisson’s ratio

σ  

=

(v P / v S  ) 2

−

2((v P / v S  ) 2

2

−

1)

Poisson’s ratio is usually around 0.5 for incompressible materials, and softer material will have values in the range 0 to 0.2 In this case, the P-wave reflection coefficient becomes larger as the angle of incidence increases. Since an increase in angle  corresponds to a greater source-receiver offset , this is called an amplitude-versus offset (AVO) anomaly. The effect of AVO effects on a synthetic CMP gather are shown with cmp_v4.m

5

6

Data examples of AVO

These effects are often present in field data and can indicate the presence of oil or gas with more reliability than normal incidence reflection data. The following example is taken from Ostrander   (1984) and shows a reflection from a known gas reservoir at a depth of 6700 feet in the Sacramento Valley, California. The amplitude increases  with offset on a number of CDP gathers.

In this example, also from Ostrander (1984), a high amplitude reflection in a sedimentary  basin in Nevada was analysed. The amplitude decreases  with offset, suggesting that the layer had a normal Poisson’s ratio. It was subsequently drilled and shown to be 160 feet of basalt.

7

C2.9.2 Images of hydrocarbon reservoirs on seismic sections

To illustrate some of the concepts used in interpretation, we will look at a number of seismic sections from a range of tectonic settings. In Labs 6 and 7 you will work though a couple of examples on your own (and realize what it takes to draw lines on a blank section!).

2.9.2.1 Extensional environment – North Sea

Brent oilfield, North Sea, Kearey Figure 4.62

Viking gas field, North Sea, Kearey 4.61. Gas reservoir is located in a faulted anticline. The most prominent reflections in a seismic section are called markers. These are correlated across the section though a comparison of their character and sequence. They are tied to lithologic units through through measuring well logs and computing computing synthetic seismograms. Vertical seismic profiles can also be used in this respect.

8

2.9.2.2 Fold and Thrust belts Wyoming, South Elk Basin Interpretation Lab 6 : Tracing reflectors across a series of faults can be more difficult in this environment, but is quite possible, as will be demonstrated in Lab 6 Rocky Mountain Foothills Turner Valley

 Anderson et al, 1989, page 165

Yan and Lines , The Leading Edge, (2001)

9

Post–stack depth migration

Pre-stack depth migration

Yan and Lines, The Leading Edge, (2001)

2.9.2.3 Reefs

Reefs are often a zone of high porosity in a carbonate layer and make good hydrocarbon reservoirs. When exposed on the surface they can form spectacular topography, such as the Permian Reef in Guadalupe Mountains National Park (Texas). However, they are less spectacular on seismic data!

Leduc reefs

● Seismic

characteristics of the reef are draping of overlying sediments and upper surface.  Note also that the base of the reef does not appear as a flat feature. It is pulled up  because of the higher velocity within the reef. 10

● Seismic

section above is from the Redwater Leduc reef 40 km northeast of Edmonton. The reservoir is Devonian age carbonates, overlain by Ireton formation shales that  provide a seal.

●Also

note the base of the upper Leduc and Duvernay appears to rise under the reef. This is an example of pull up. The base is flat and appears at earlier TWTT since the velocities in the reef are higher than at the same depth east of the reef. This is also  Andersen et al., 1989, Chapter 1.) illustrated in the synthetic example on the left ( Andersen

11

● 

Again, well logs are used to generate synthetic seismograms, which allow the interpretation of the seismic section. This will show which lithologic contrasts will dominate the seismic response. Allows the structure mapped in a well to be extrapolated away from the well in places where seismic data are available.

Louisiana and East Texas reefs

Telford 4.104, Telford 4.105

2.9.2.4 Salt related hydrocarbon reservoirs reservoirs

Salt can produce significant hydrocarbon traps when sedimentary units are deformed in salt tectonics. Some of the earliest hydrocarbon deposits discovered with geophysics were associated with salt-domes, and were detected through gravity surveys with torsion balances and pendulums (Geophysics 224, Section B)

12

Salt creeps under mechanical loading and salt sheets can flow horizontally. Salt diapers can rise in a column to the surface. The figure above shows 4 second seismic reflection data that reveals the evolution of a salt sheet, progressing from a salt swell, to a diapir that  penetrates the sediments, and finally to the surface (North Sea, Telford figure 4.102)

Hydrocarbon reservoirs may be located in the sedimentary sequence beneath a salt layer. However, the high velocity of a salt body makes it difficult to image these reservoirs.

The examples above show the Enchilada, Hickory and Tanzanite sub-salt discoveries in the Gulf of Mexico, and are taken from an advert by Diamond Geophysical Service Corporation ( The Leading Edge, November 1999) Brown denotes the salt body and yellow the reservoir. Imaging techniques such as pre-stack depth migration have made a  big improvement in this type of exploration problem, as described in C2.8 for the Sigsbee dataset.

3-D seismic reflection surveys

In a 3-D seismic survey data is recorded on a grid of surface lines. In land, the shot lines and recording lines are generally orthogonal. ( Kearey 4-34). This produces a grid of CDP’s on a subsurface reflector.

In marine seismic surveys, an array of airguns and multiple streamers give a swarth of 3D coverage (Kearey 4-36). Once 3-D data has been collected, it may displayed as a seismic data volume (seismic cube). A vertical slice corresponds to the depth sections we have discussed.

13

Three-dimensional seismic data volume (seismic data cube) collected over a salt dome in the Gulf of Mexico. Kearey Plate 4.1. Colour corresponds to amplitude of reflection.

A time slice (seiscrop section) from the above data volume at 3.7 seconds. At each CDP the amplitude of the trace is  plotted in colour. Note that changes red black correspond to positive-negative arrivals (layers). Steeply dipping layers  produce closely spaced +/- changes.

Time-structure Time-structure maps – maps depth to reflector in time ( e.g. Kearey 4.47 ) Structural contour maps - maps depth to reflector in depth  Isochron maps – shows variation of layer thickness as function of time thickness in metres  Isopach maps maps - layer thickness

Other interpretational techniques and concepts Seismic stratigraphy – relate the seismically imaged sedimentary sequence to depositional history. Locate unconformities etc. Seismic facies analysis. Use of the reflection geometry and character to determine the sedimentary facies. This can reveal the environment in which sedimentary rocks were deposited. Can help identify potential source rocks, seals and reservoir rocks.

14

Sometimes 3-D surveys can show paleo river channels in time slices. In this location a particular reflector is deeper over the river channel and at a given time, this will give a different reflection amplitude on CDP’s that are within/outside the channel. Kearey Plate 4.4

C2.9.3 Recent developments developments in hydrocarbon exploration exploration with seismic reflection

2.9.3.1 Pre-stack depth migration . ●

 Now possible with the increased computing power that is available for processing seismic data. With steeply dipping structures, this technique is much better at imaging than post-stack depth migration schemes. The Leading Edge, 21, Special section PSDM 2, The Sequel, December 2002.

2.9.3.2 Marine exploration ●

Concerns about the possible effects of air gun noise on marine mammals has intensified in recent years. Details described by Gausland , 2000.

2.9.3.3 Multiple attenuation.

Multiples can obscure later arrivals in reflection data. Various processing techniques can  be used to remove them. Matson et al., (1999) give a comparison of various methods.

 - Multiples in data from the Canarvon basin (Australia) from  Matson et al., (2001).  Left   – section after the most convincing multiple suppression technique was applied.  Right  – Vertical scale is 0 to 1.2 seconds

15

Pre-stack depth migrated images  before and after multiple removal. Vertical scale is 1200 to 6800 m. Subsalt geology is definitely visible after multiple removal. From Guitton (1999).

● what

will be the repeat times of multiples in 1 km of water? In 2 km? vw= 1500 m/s

● S-wave

data can have less problems with multiples than P-wave data? Why?

2.9.3.4 Amplitude versus offset ●

State-of-the-art :The Leading Edge, 19, Special section AVO, November 2000.

●

AVO doesn’t always work! Some cautionary tales from Allen et al, (1993)

2.9.3.5 Shear waves and seafloor recording ●

3-component (3-C) recording allow S-waves to be recorded included in the analysis. This requires seafloor  seismic recording in marine surveys since S-waves cannot travel through the water column. The converted phase is sometimes called a C-wave.

● Use

either ocean bottom seismometer (OBS) or ocean bottom cables (OBC). These can be laid in a fixed location on the seafloor, or dragged by the survey ship.

  S-waves are generated by P-to-S conversion at the reflector. This requires nonnormal incidence, so large source-receiver offsets are used.

●

16

Imaging through gas clouds

Seismic reflection image of an anticline in the Tommeliten Field, Norwegian sector of of the North Sea. Gas escaping from the reservoir forms a gas chimney, which attenuates the P-waves. This occurs because P-waves are sensitive to the compressibility of gas in pore spaces. In contrast, S-waves do not compress the rock as they travel through a rock unit. From Caldwell (1999)

The same structure imaged with the arrival of S-waves produced by P-to-S conversion at depth. These waves are significantly less attenuated by the gas cloud and the top of the anticline is well imaged.

Some reservoirs have no P-wave expression

Towed streamer data (P-waves)

Ocean bottom cable (S-waves)

In the Alba Field (North Sea) the reservoir sands are overlain by a shale cap rock. The change in P-wave velocity is small at the top of the reservoir, and it is not imaged on the P-wave section. However there is an increase in S-wave velocity that produces a strong reflection (Caldwell, 1999).

17

Better illumination of high angle interfaces ● 

Seafloor recording of P-waves and S-waves can also give wider angle data coverage. This can improve illumination of targets. Example from Mississippi canyon, Gulf of Mexico, PGS advert, The Leading Edge, 20, 28, January 2001

Towed streamer

Seafloor cable

The Leading Edge, 18 , Special section : Multi-component offshore , November 1999 1999 The Leading Edge, 20 , Special section : Advances in shear wave technology , September 2001.

2.9.3.6 Detecting time variations  in hydrocarbon reservoirs during production ● Repeat

seismic survey at same location to determine time variations during production (time lapse or 4-D seismics)

● 

Need very careful processing to make an objective comparison between surveys. Ideally use exactly the same source and geophones. In some cases, the geophones are left in place between surveys. Cross-equalization is a processing technique that looks at events in trace that are not associated with the reservoir and seeks to equalize them over time.

Example from steam flood of a heavy oil reservoir at East Senlac SK. The survey was repeated in 1990, 1997 and 1998. Unfortunately different seismic sources and geophones were used in each survey, so an objective comparison is complicated (Li et al., 2001). These results suggest that the negative  polarity event at the top of the reservoir  becomes stronger over time as steam heats the reservoir (maybe????)

18

●  What

other factors than the reservoir could change seismic data recorded at the same location in different surveys?

●

Differencing of seismic sections. Kearey 4.41 shows time sections before and after gas gas was injected into a reservoir for storage.

●  Kearey

4.42. Steam injection at a well in Duri oilfield, Indonesia. Base of reservoir initially shows pull-up and then push-down (sag). What does this tell us about the velocity in the reservoir during the injection?

References

Allen, J., C. Peddy and T. Fasnacht, Some AVO failures and what (we think) we have learned, The Leading Edge, 163-167, March 1993. Andersen, N.L. et al, Geophysical Atlas of Western Canadian Hydrocarbon Pools, Canadian Society of Exploration Geophysicists, 1989. Caldwell, J., Marine multi-component seismology, The Leading Edge, 18, 1274-1282, 1999. Gausland, I., Impact of seismic surveys on marine life, The Leading Edge, 19, 903, 2000. Guitton, A., Multiple elimination using a pattern recognition technique, The Leading  Edge, 18, 92-99, 1999. Li, G. et al, Effective processing of non-repeatable 4-D seismic data to monitor heavy oil SAGD steam flood at East Senlac, Saskatchewan, Canada, The Leading Edge, 20, 5459, 2001. Matson, K., D. Paschal, A. Weglein, A comparison of three multiple-attenuation methods applied to a hard water-bottom data set, The Leading Edge, 18, 120-126, 1999. Macleod, M.K., et al, The Alba Field ocean bottom cable survey: Impact on development, The Leading Edge, 18, 1306-1312, November 1999. Ostrander, W.J., Plane-wave reflection coefficients for gas sands at non-normal angles of incidence, Geophysics, 49, 1637-1648, 1984. Telford, W.M., L.P Geldart, R.E. Sheriff, Applied Geophysics, Cambridge University Press, 1990.

Geophysics 224 MJU 2006 

19

C2.10 Deep seismic reflection studies

In addition to their use their  use in hydrocarbon exploration, seismic reflection studies can also be used to image the entire crust and upper mantle. A number of government and university research groups have made a number of deep seismic studies in recent years and improved our knowledge of the structure and evolution of the Earth. 2.10.1 BIRPS (British Institutions Reflection Profiling Syndicate)

The first example is from the British Institutions Reflection Profiling Syndicate (BIRPS) and is shown in Kearey Figure 4.63. BIRPS collected most of their data with marine seismic surveys around Great Britain where the continental crust is present on the continental shelf. This reflection section is just like the sections collected for hydrocarbon exploration, except the recording time is longer (16 seconds). Collection of these data require a strong seismic source, typically a very large airgun array. Often more energy was used by the compressors than the ships engines. A very long long seismic streamer was used to detect the normal move out in deep reflections. Some BIRPS surveys used to 2 ships to create a longer offset.

This BIRP clearly shows the Moho at about 8 seconds and a major thrust fault that apparently extends into the upper mantle. A common feature of many BIRPS profiles is significant seismic reflectivity in the lower crust. The origin of this feature is still not defined, but lenses of fluids, or mafic underplating are two possible explanations.

1

2.10.2 COCORP (Consortium for Continental Reflection Profiling)

The Consortium for Continental Reflection Profiling (COCORP) used land seismic data to study deep crustal structure in the United States.

The figure above shows a seismic reflection image of the Wind River Thrust in Wyoming. This is a major feature associated with the Laramide orogeny ( Brewer et al, 1980). Typical Typical acquisition parameters for COCORP data were VIBROSEIS with a 96fold CMP stack, recorder spacing of 100 m and spread length of 9.6 km. The VIBROSEIS signal was 8-32 Hz and at each location was stacked up to 80 times.

2.10.3 INDEPTH (International Deep profiling of Tibet and the Himalaya)

By the early 1990’s COCORP had examined the crustal structure of many regions of the United States. This included old orogenic belts such as the study above, and active features such as the San Andreas Fault and Basin and Range province. However not all tectonic processes are currently active in the United States and they began to look overseas to places where certain tectonic processes are active today. Interest in continent-continent collisions led to the INDEPTH project in Tibet. This used a Chinese seismic crew to collect seismic reflection and refraction data on a profile that crossed the Himalaya at 27.5º N and extended north onto the Tibetan Plateau. The  primary energy source was explosives.

2

These seismic data imaged the Indian Plate being thrust beneath the Himalaya. The Main Himalayan Thrust (MHT) is the top of the Indian Plate. Note also that the Moho in the Indian Plate can be seen at a two-way travel time of 25 seconds. This observation showed that the 80 km thick crust was made by stacking Asian and Indian crust on top of each  Nelson et al, 1996). other ( Nelson Another key observation from INDEPTH was that bright spots  were observed beneath the Tibetan Plateau extending north from 29.75º N. These have a high amplitude and negative polarity, and may represent a thin layer of fluids.

3

Panel on the below shows how the polarity of the reflection (YBS) is determined by comparison with the first break (FB) of the direct P-wave. Since explosives were used as  Brown et al., a source, it is known known that the signal has a compression as the first arrival ( Brown 1996).

Combined with other seismic and petrologic data, many researchers believe they represent a zone of partial melt, overlain by aqueous fluids. Additional evidence for widespread melting comes from magnetotelluric data that show that low electrical resistivity is coincident with the seismic reflections (Unsworth et al., 2005).

4

2.10.4 LITHOPROBE

Lithoprobe was a federally funded geoscience project that used combined geology and geophysics to study the structure and tectonic evolution of a number of regions of Canada  between 1980 and 2000. Deep seismic reflection reflec tion was an integral part of several transects. This included marine seismic data collection on the west coast and in the Great Lakes. Onshore seismic data was collected in other transects, including the Alberta Basement Transect.

The map shows the location of seismic reflection lines in Alberta (black lines). A typical depth section is shown below that extends north-south across the Vulcan structure. The Moho is visible on the right side at a depth of 35-40 km. More details of interpretation are described by Bouzidi et al., (2002).

5

References

Bouzidi, Y., D. R. Schmitt, R. Burwash, E. Kanasewich, Depth migration of deep seismic reflection profiles: crustal thickness variations in Alberta, Can. J. Earth Sciences, 39, 331-350, 2002. J. A. Brewer, S. B. Smithson, J. E. Oliver, S. Kaufman and L. D. Brown, The Laramide orogeny: Evidence from COCORP deep crustal seismic profiles in the wind river mountains, Wyoming, Tectonophysics, 62, 165-189,1980. Brown, L. D., W. Zhao, K. D. Nelson, M. Hauck, D. Alsdorf, A. Ross, M.Cogan, M. Clark, X. Liu, J. Che, Bright spots, structure and magmatism in southern Tibet from INDEPTH seismic reflection profiling, Science, 274, 1688-1690, 1996.  Nelson, K.D., Zhao, W., Brown, L.D., Kuo, J., Che, J., Liu, X., Klemperer, S.L., Makovsky, Y., Meissner, R., Mechie, J., Kind, R., Wenzel, F., Ni, J., Nablek, J., Leshou, C., Tan, H., Wei, W., Jones, A.G., Booker, J.R., Unsworth, M., Kidd, W.S.F., Hauck, M., Alsdorf, D., Ross, A., Cogan, M., Wu, C. Sandvol, E., Edwards, M., Partially molten Middle Crust Beneath Southern Tibet : Synthesis of Project INDEPTH results , Science, 274, 1684-1686, 1996. Unsworth, M.J. A.G. Jones, W. Wei, G Marquis, S. Gokarn, J. Spratt, Crustal rheology of the Himalaya and Southern Tibet inferred from magnetotelluric data, Nature, 438, 7881, doi:10.1038/nature04154, doi:10.1038/nature04154, 2005.

6

C2.11 Vertical seismic profiling

Kearey Figure 4.44

1

Applications of vertical seismic profiling 1. Time-depth calibration

What depth are reflections coming from? Primary and multiple identification Does an interface actually produce a reflection?

2. Higher resolution imaging

Surface 3-D

VSP

2

C3.1 Seismic refraction refraction – single horizontal horizontal interface interface

The P-wave is refracted at the interface between the two layers. Since v1 > v2 the wave is refracted towards   the horizontal. As the angle of incidence is increased, the geometry results in a head wave  travelling horizontally  in layer 2. From Snell’s Snell’s Law we can write:

sin θ  c v1 Thus sin θ c

=

v1 v2

=

sin 90

o

v2

and from geometry we can show that AB = CD =  z / cos θ c

Can also show that BC = x =  x - 2z tan θ c Total travel time

t = tAB t  = t  = t  = t  = t  =

t  = t  =

+ tBC  z 

v1 cos θ c  x v2  x v2  x v2  x v2  x v2  x v2

+

+

+

+

+

+

+

+ tCD

( x − 2 z tan θ c )

2 z  v1 cos θ c

v2 −

+

 z  v1 cos θ c

2 z tan θ c v2

2 zv2 − 2 zv1 sin θ c v1v2 cos θ c 2 z − 2 z sin 2 θ c v1v2 cos θ c 2 z cos θ c v1 2 z  v22 − v12 v1v2 constant

1

•

The travel time curve for the refracted wave is a straight line with slope = 1 / v2

•

The refracted arrival is only observed when x when  x > xcrit = 2z tan θ c

•

The refracted wave is the first arrival at values of  x greater than the cross over distance ( x  xcross)

•

When x When x = xcrit  the refracted and reflected waves are the same

•

v2 can be calculated from the slope of the refracted wave on the t-x plot t-x plot

•

The depth of the interface ( z ) can be found by extrapolating the travel time of the refracted wave to x to x =  = 0 where the travel time is

t i

Rearranging gives

=

 z  =

2 z  v 22

−

v12

v1v 2 v1v 2 t i 2 v 22

−

v12

2

C3.2 Seismic refraction refraction – multiple horizontal horizontal layers

Direct wave

t  =

 x v1

 Data analysis ● compute v1 from slope of direct arrival

First refraction

Critical refraction occurs as the wave travels from 1>2 giving v  sin θ c = sin θ 12 = 1 v2

t  =

 x v2

+

2 z 1 cosθ 12 v1

=

 x v2

t 

+ 1

 Data analysis ● Compute v2 from the slope of the refracted wave  Compute z 1 from the intercept time (t   (t 1) when x when x =0 and v1 and v2 are already known ● Compute z  2 z  cosθ 12 t 1 = 1 v1

Second refraction

Critical refraction occurs as the wave travels from 2 > 3 giving sin θ 23 Applying Snells Law at the interface between 1 and 2 gives sin θ 13

=

=

v1 v2

v2 v3 sin θ 23

3

Thus sin θ 13

=

v1 v2

sin θ 23

=

v1 v3

Total travel time for the second refraction  x 2 z 1 cosθ 13 2 z 2 cosθ 23  x + + = t  = v3 v1 v2 v3

t 

t 

+ 1 + 2

 Data analysis ● compute v3 from the slope of the second refracted wave ● At x  At x =0, =0, intercept time = t 1+t 2 ● Since t 1 is known we can compute t 2 ●  From t 2

=

2 z 2 cosθ 23 v2

 we can then calculate z  calculate z 2

Example of field seismogram

4

Multiple layers with velocity increasing with depth

The results for two layers can be generalized for an N  an  N -layer -layer model and written as a series. See the text book for details. Consider the case where there is a uniform increase in velocity with depth. What will be the form of the ray paths and travel time curves?

What could cause velocity to increase uniformly with depth? For each ray, the quantity called the ray parameter , p,  p, is constant. constant. p =

sin θ  v

Low-velocity layers

● for a head wave to propagate, an increase in velocity  from one layer to the next is needed. If a decrease in velocity occurs, there will be no head wave and refraction will fail to detect the layer.  Example :  Example : If a 3-layer model has a LVZ in the second layer, then interpretation in terms of 2-layers will give a wrong answer e.g. soil-peat-bedrock. The peat is a low velocity layer. Hidden layers

Even if a layer has an increase in velocity, then it possible for the head wave on the upper surface to never be the first arrival . Again if a 3-layer velocity model has a hidden second layer, then interpretation in terms of 2-layers will give a wrong answer

5

6

7

C3.3 Seismic refraction refraction – dipping dipping interface interface Direct wave

t  =

 x v1

● same as in horizontal case considered in C3.2 ● compute v1 from slope of direct arrival

Refraction from the dipping interface

If ray is travelling down dip, then the upward leg to the geophone will increase in length as offset ( x)  x) increases. Thus the first refraction will arrive  progressively later, compared to the flat interface in C3.1. This if effectively the same as the signal travelling more slowly and will increase the slope  of the refracted wave on the travel time curve ( x-t   x-t  plot)  plot)

● Can show that the travel time, t  for 

=

 x sin(θ 12

+

γ 1 )

v1

+

2 z  perp cosθ 12 v1

● Can you confirm that this result is correct when γ1 = 0? ● Apparent velocity down dip

=

v 2 d 

=

v1 sin(θ 12

+

γ 1 )

< v2

● Note that the apparent velocity down dip is slower  than v2 which results in a steeper slope on the travel time curve.

● Example of layer with γ1 = 5 v1 = 1000 m/s and v2 = 3000 m/s is shown on page 9. ˚

What value of v 2d is predicted for these values?

8

9

● Note that several velocity models can be found found that all fit the data. v1  is uniquely determined, but there is a trade-off between γ1 and v2 (non-uniqueness again)  Model 1  Model 2  Model 3

γ1= 0 v1 = 1000 m/s v2 = 2400 m/s γ1= 5 v1 = 1000 m/s v2 = 3000 m/s γ1=10 v1 = 1000 m/s v2 = 3900 m/s ˚

˚

˚

10

●  On the basis of a single shot recorded in one direction, we cannot determine which model is the correct one. ● Suppose we fire a second shot at the far end of the geophone array, then these waves will travel up dip. As the refraction travels to larger offsets ( x)  x) the final leg in the upper layer will become shorter, and the refraction will arrive earlier. This effectively increases the apparent velocity  and reduces the slope  of the travel time curve.

t rev

=

 x sin(θ 12

−

γ 1 )

v1

● Apparent velocity up dip

+

=

' 2 z  perp cosθ 12

v1 v 2u

=

v1 sin(θ 12

−

γ 1 )

For with γ1 = 5 v1 =1000 m/s and v2 = 3000 m/s, what value of v 2u is predicted? ˚

●  Note Note that the travel times tAD  (forward direction) and t DA  (reverse direction) are the same. This is a phenomena known as reciprocity. Why is this true? ● When forward and reverse profiles are recorded, the correct model (#2) can be seen to  be the only one that fits the travel time data (attached handout) Computing v2 and dip of the interface

By rearranging expressions for v 2u and v2d can show that

θ 12

=

1 2

[sin 1 ( −

v1

) + sin 1 ( −

v 2 d 

v1 v 2u

)]

and

γ 1

=

1 2

[sin 1 ( −

v1 v 2 d 

) − sin 1 ( −

v1 v 2u

)]

11

12

C3.3 Seismic refraction refraction – dipping interface interface - Example Example

Look at the travel time curve. Which way does the interface dip? You will need this information to decide which direction is up and which is down! down! Forward profile x(m) t _di _di r ( ms) 0. 000 000 0. 000 000 50. 50. 000 000 33. 33. 000 000 100. 100. 000 67. 000 150. 150. 000 100. 100. 000 200. 200. 000 133. 133. 000 250. 250. 000 167. 167. 000 300. 300. 000 200. 200. 000 350. 350. 000 233. 233. 000 400. 400. 000 267. 267. 000 450. 450. 000 300. 300. 000 500. 500. 000 333. 333. 000 550. 550. 000 367. 367. 000 600. 600. 000 400. 400. 000 650. 650. 000 433. 433. 000 700. 700. 000 467. 467. 000 750. 750. 000 500. 500. 000 800. 800. 000 533. 533. 000 850. 850. 000 567. 567. 000 900. 900. 000 600. 600. 000 950. 950. 000 633. 633. 000 1000. 1000. 000 667. 667. 000

t _r ef ( ms) 110. 110. 410 410 133. 133. 927 927 157. 157. 444 444 180. 180. 960 960 204. 204. 477 477 227. 227. 994 994 251. 251. 510 510 275. 275. 027 027 298. 298. 544 544 322. 322. 060 060 345. 345. 577 577 369. 369. 094 094 392. 392. 610 610 416. 416. 127 127 439. 439. 643 643 463. 463. 160 160 486. 486. 677 677 510. 510. 193 193 533. 533. 710 710

Reverse profile x(m) t _di _di r ( ms) 0. 000 000 667. 667. 000 000 50. 000 633. 633. 000 000 100. 100. 000 600. 600. 000 000 150. 150. 000 567. 567. 000 000 200. 200. 000 533. 533. 000 000 250. 250. 000 500. 500. 000 000 300. 300. 000 467. 467. 000 000 350. 350. 000 433. 433. 000 000 400. 400. 000 400. 400. 000 000 450. 450. 000 367. 367. 000 000 500. 500. 000 333. 333. 000 000 550. 550. 000 300. 300. 000 000 600. 600. 000 267. 267. 000 000 650. 650. 000 233. 233. 000 000 700. 700. 000 200. 200. 000 000 750. 750. 000 167. 167. 000 000 800. 800. 000 133. 133. 000 000 850. 850. 000 100. 100. 000 000 900. 900. 000 67. 000 000 950. 950. 000 33. 000 000 1000. 1000. 000 0. 000 000

t _r ef ( ms) 533. 533. 710 710 517. 517. 616 501. 501. 522 485. 485. 428 469. 469. 334 453. 453. 240 437. 437. 146 421. 421. 051 404. 404. 957 388. 388. 863 372. 372. 769 356. 356. 675 340. 340. 581 324. 324. 487 308. 308. 393 292. 292. 299 276. 276. 205 260. 260. 111 244. 244. 017

13

From the travel times compute the following. v1

= __________ m/s

v2d

= __________ m/s

v2u

= __________ m/s

θ 12 12

= __________ degrees

v2 = v1 / sin θ 12 12

= __________ m/s

γ1

= __________degrees __________degrees

Use the modelling program refract_v3.m to confirm the answers.

14

C3.4 Seismic refraction – non planar interfaces

C3.4.1 Basics and concept of delay time ●

the delay time at the shot is the extra time needed for the wave to travel AB, compared to the time to travel CB.

●

delay time at detector is the extra time needed for wave to travel DE, compared to DF

●

can see this from the travel time curve

●

for horizontal interface, δts = δtd = z  = z cos θ/v1

●

with just a single shot, we cannot separately determine δts and δtd

●

with a reversed profile, we can individually determine δts and δtd

t = x/v2 + δts+ δtd

●  for a non horizontal interface, with many geophones δtd  contains information about how the depth depth of the interface (z) varies with position (x’) ●

two common methods of finding δtd are the plus-minus method and the generalized reciprocal method. More details shortly.

15

C3.4.2 Plus-minus interpret i nterpretation ation method

Consider the model with two layers and an undulating interface. The refraction profile is reversed with two shots (S1 and S2) fired into each detector (D). Consider the following three travel times:

(a) The reciprocal time is the time from S 1 to S 2

t S 1S 2

=

(b) Forward shot into the detector 

t S 1 D

=

(c) Reverse shot into the detector

t S 2 D

=

l  v2  x v2

+

+

δ S 1

δ S 1

(l − x) v2

+

+

+

δ S 2

=

t S 2 S 1

δ  D

δ S 2

+

δ  D

Our goal is to find v2 and the delay time at the detector, δ D. From the delay time, δ D , we can find the depth of the interface.

16

Minus term to estimate velocity (v 2)

(b)-(c) will eliminate δ D

where C  is  is a constant. A plot of t S  D 1

t 

− S   D 2

t S 1 D

− S   D = 2

t 

t S 1 D

− S   D = 2

t 

( 2 x − l ) v2 2 x v2

+

+

δ S 1

−

δ S 2

C 

 versus 2 x will  x will give a line with slope = 1/v2

Plus term to estimate delay time at the detector

(b)+(c) gives

t S 1 D

Using the result (a) we get

t S 1 D

Re-arranging to get an equation for δD

δ  D

t 

+ S   D = 2

t S 2 D

+

=

1 2

=

(t S 1 D

l  v2

+

t S 1S 2

+

δ S 1

+

t S 2 D

+

δ S 2

+

2δ  D

2δ  D

t 

− S  S  1 2

)

This process is then repeated for all detectors in the profile

17

Example 1 : Synthetic refraction data were generated for a dipping flat interface with the MATLAB script refract_v3_data_plus_minus.m

The plus-minus technique plus_minus_ex1.m

was

then

applied

using

the

MATLAB

Original model 

v1=1500 m/s

v2=2500 m/s

 Interpretation with plus-minus method 

v1=1502 m/s

v2=2525 m/s

program

The interface dips to the right at 8 , and is correctly imaged in this example. ˚

18

Example 2: Field data from Kearey and Brooks, p.123

19

C3.4.3 Generalized Generalized reciprocal method

The plus-minus method assumes a linear interface between points where the ray leaves the interface. A more powerful technique if the Generalized reciprocal method in which  pairs of rays are chosen that leave the interface at the same location. More details can be found in Kearey in Kearey page 109. 109.

C3.5 Shallow applications of seismic refraction

C3.5.1 Depth of bedrock

•

velocity of bedrock greater than unconsolidated layer

•

in this example, (Kearey 5.24) a shot point was located every 30 m

•

depth to bedrock increases with x with  x

20

Example from Northern Alberta

Seismic refraction was used to determine depth to bedrock at the location where a  pipeline was planned to cross a creek. DC resistivity data were also collected at this location. Data courtesy of Greg Kovacs and Komex

 Note that the direct wave is the only the first arrival at the first 2 geophones. This is  because of a very high velocity contrast between the upper and a nd lower layers. The model  below was derived from the seismic data using the general reciprocal method.

21

3.5.2 Locating water table

3.5.3 Rippability

22

C3.6 Crustal scale seismic refraction surveys

3.6.1 Discovery of the crust-mantle boundary (the Moho)

In 1910 Andrija Mohorovicic published a paper describing the time taken for earthquake waves to travel across the Balkans. He noticed that P-waves apparently travelled at two different velocities. He inferred that this required a crustal layer, underlain by a higher velocity layer. This boundary is now called the Mohorovicic discontinuity or Moho.

●What crustal velocity and crustal thickness is implied by the attached figure? ●How do the travel time curves differ from those we have considered with a hammer or explosive source? ●Why can the S-waves be detected, since they are not first arrivals?  P  = P  g  ● In modern seismic terminology  P = P n and  P =

23

24

● the direct wave is referred to as Pg and often travels as a parabola. Why? ● the head wave that travels in the upper mantle is called Pn ● reflection from the Moho is called P mP ● reduced travel time  is sometimes  plotted on the vertical axis. t' = t – x/v red where vred   is the reduction velocity . This has the effect of making arrivals with v=vred   plot horizontally  on a t-x  plot.

● in the figure on the left, the crustal Pwave velocity was used as the reduction velocity.

3.6.2 Seismic refraction at sea

●  signals can be detected either by seafloor recording, using a long streamer, two-ship survey or a sonobuoy. ●  refractions (head waves) occur at seafloor and at deeper interfaces. ●  sonobuoy sends data back via radio telemetry. Modern sonobouys are often disposable.

25

●  Can you locate the direct water wave on this section? (P-wave velocity is 1500 m/s) ●  Locate a refraction and measure the velocity. ● Where would you expect to see normal incidence reflections?

3.6.3 Fan shooting

Technique first used in the 1920’s in the search for salt domes. The higher velocity of the salt causes earlier arrivals for signals that travel though the salt.

Kearey and Brooks

Eve and Keys, Applied Keys, Applied Geophysics, Geophysics, 1928

26

3.6.4 Tectonic studies of the continental lithosphere with refraction Deep Probe experiment

●

used explosive shots up to 2400 kg with seismic recorders deployed on a profile from 60  N to 43  N ˚

˚

● ray tracing used to model the data. Measures the variation in Moho depth and crustal structure. Figure above shows ray tracing for SAREX shot 5. Note that with a reduction velocity of 8 km/s, P n plots as a horizontal line, while the slower Pg has a positive slope. Gorman, A.R. et al , Deep probe: imaging the roots of western North America, Canadian  Journal of Earth Sciences, Sciences, 39, 375-398, 2002. Clowes, R. et al , Crustal velocity structure from SAREX, the Southern Alberta seismic experiment, Canadian Journal of Earth Sciences, Sciences , 39, 351-373, 2002.

Geophysics 224 MJU 2006 

27

Geophysics 224 C. SEISMIC EXPLORATION METHODS- Outline Basic principles of seismology 1.1

General introduction to wave phenomena Frequency, wavelength, period. Phase. Angular frequency. Ray paths and wave fronts. Huyghens principle.

1.2

Wave propagation in Earth Body waves - P-wave, S-wave, Surface waves – Rayleigh and Love waves Particle motion in movie. Stress and strain. Elastic modulus, shear modulus. Poisson’s ratio. v

k  =

 ρ 

1.3

Typical seismic velocities for Earth materials Typical range of values. P-wave velocity in porous rocks. Time-average equation. Velocity as a function of depth Rippability

1.4

Reflection coefficients at normal incidence Impedance, Z = density x velocity

For amplitude

 R

T 

=

=

v 2 ρ 2

−

v1 ρ 1

v 2 ρ 2

+

v1 ρ 1

2v1 ρ 1 v 2 ρ 2

+

v1 ρ 1

Energy partitioning 1.5

Reflection and refraction at oblique incidence. Snell’s Law. Critical angle Head waves

1.6

Diffraction Migration is the solution.

1.7

Factors that cause the amplitude to change as wave propagates. Geometrical spreading. Er 2 = constant. Ar = constant. Attenuation of seismic signals and exponential decay. Frequency dependence causes dispersion. Decibels.

1.8

Seismic energy sources Air guns, explosives Vibroseis, weight drop, hammer, Betsy gun. Stacking of signals.

1.9

Seismic detectors. Geophones and hydrophones

2. Seismic reflection 2.1

Travel time curve for single horizontal interface.  Normal moveout (NMO) to separate depth and velocity.

2.2

Travel time curves for multiple horizontal layers. Dix equation. Multiples.

2.3

Travel time curves for a dipping reflector.

2.4

Displaying real seismic data. Traces and seismic sections. Reduction velocity.

2.5

Common depth-point profiling. Stacking to improve signal-to-noise ratio.

2.6

Resolution. Vertical limitation from wavelength. Horizontal limitation from size of Fresnel zones.

2.7

Static corrections and sources of noise. Filtering strategies to reduce noise. Convolution and deconvolution.

2.8

Migration to recover true reflector geometry and depths.

2.9

Applications-I. Hydrocarbon exploration. Structural analysis. 3-d acquisition. Time slices. Seismic stratigraphy.

2.10 Applications-II. Lithospheric studies. BIRPS, Lithoprobe and INDEPTH data. 2.11 Vertical seismic profiling.

3. Seismic refraction 3.1

Travel time curves for a two layer Earth.

3.2

Travel time curves for a three layer Earth. Hidden layers (thin or low velocity)

3.3

Travel time curve for a dipping layer. Reversed refraction profile.

3.4

Applications in shallow geotechnical studies

3.5

Applications in lithospheric studies. Discovery of the Moho. MJU 2003

Geophysics 224 D1: Basics of geomagnetism

D1.1 Definitions H = magnetic field strength B = magnetic flux density μ =

magnetic permeability

B = μ H

In the absence of magnetic materials, μ = μ0 = 4π x 10-7 H/m

D1.2 Magnetic monopoles

Despite extensive searches, magnetic monopoles do not appear to exist in isolation. They always occur in pairs of positive and negative monopoles that form dipoles (see D1.4).

Breaking a bar magnet (a dipole) in half does not generate separate monopoles, rather two new dipoles.

http://www.oberlin.edu/physics/catalog/de http://www.oberlin.edu/physics/catalog/demonstrations/em/mag monstrations/em/magneticmonopole.html neticmonopole.html

Some theories in particle physics predict that monopoles may be observed in high energy collisions between subatomic particles: http://www.aip.org/png/html/monopole.htm

D1.3 Magnetic field lines

Many animals have the ability to detect the direction of the Earth’s magnetic field and use it for navigation. Note the recent study of pigeons by Mora et al., (2004) listed below. However most  people cannot detect the magnetic field fie ld of the Earth and other ways are needed to visualize the magnetic field.

Mora, C.V., Davison, M., Wild, J.M. and Walker, M.M. Magnetoreception and its trigeminal mediation in the homing pigeon,  Nature, 432:508-511, 2004. Magnetic field lines represent the direction in which a magnetic pole would move. The force, F, on a pole of strength m is defined as H = F/m

D1.4 Magnetic dipoles

Consider a magnetic dipole, with poles m+ and m-, separated by a distance l. The magnetic dipole moment is defined as P = ml The total magnetic field is the vector sum of m+ and m-

D1.5 Diamagnetism and paramagnetism

On the atomic level, magnetic behaviour is due to atoms behaving as small magnetic dipoles. If a uniform magnetic field (H) is applied, there are two possible responses. Diamagnetic behaviour

The when a magnetic field ( H) is applied, the atoms develop a dipole moment that opposes the applied magnetic field. The resulting magnetic moment (M) is related to the applied magnetic field (H) as M=k H where k  is   is defined as the The magnetic susceptibility. magnetization (M) is said to be induced magnetization, since it will disappear when the field H is removed. For a diamagnetic material, k  is  is small and negative. Diamagnetic materials commonly found in the Earth include salt, quartz and feldspar.

Paramagnetic behaviour

The atoms align parallel to H and increase the magnetic field. B

= μ (H+M) = μ (H+k H) = μ (1+k) H = μ μr  H

where μr   is defined as the relative permeability. For paramagnetic materials  is positive. k  is Paramagnetic elements include iron, nickel and cobalt. Again M  is described as induced magnetization   since it will vanish when H is removed.

D1.6 Ferromagnetism and magnetic domains

In certain minerals the  paramagnetic behaviour is especially strong. If a few atoms  become aligned with an applied magnetic field, then the field in the material increases and more atoms  become aligned. Through positive feedback a whole region of the material can become magnetized, and is termed a domain. This spontaneous magnetization does not disappear when the applied field is removed and it is termed remnant magnetization. Domain theory accounts for observations of magnetization in real minerals, and paramagnetic  behaviour can be subdivided into ferromagnetic  and ferrimagnetic  behaviour. The size and orientation of the domains is determined by lattice structure.

Haematite: Magnetite:

Magnetic domains in garnet

Magnetic stripe domains in Co-Zr-O2

anti-ferromagnetic  ferrimagnetic

A very simplistic computer simulation  of ferromagnetic behaviour is illustrated in the MATLAB script paramag.m  and shows a set of atoms that have a magnetic moment and which are free to rotate. When an external magnetic field is applied, the atoms begin to line up. The magnetic field of each atom can influence its neighbours and complex patterns can develop. Example 1 shows a small grid (11 x 19) of dipoles (atoms). Note that the internal magnetic field strength becomes strong in regions where ordering occurs. Example 2 shows a 300 x 300 grid.

D1.7 Curie temperature t emperature

As temperature increases, thermal vibration energy begins to breakdown the ordering of a ferromagnetic material. Above the Curie temperature, spontaneous magnetization ceases.

Typical Curie temperatures: Pure iron

1043 K

Fe2O3

893 K

For a lab demonstration see http://www.geol.binghamton.edu/faculty/bark http://www.geol.bingha mton.edu/faculty/barker/demos/demo er/demos/demo13.html 13.html There are two important consequences of the Curie temperature (1) Since temperature in the Earth increases with depth, there exists a depth below which materials cannot behave as ferromagnetic. Thus only rocks at shallow depths in the Earth can exhibit remnant magnetization. (2) As a rock cools from above the Curie temperature, it will acquire a magnetic field that records the strength and direction of Earth’s field at that time. By using radioactive

dating to find the age of the rock, this gives us a powerful tool to determine how the Earth’s magnetic field has varied over time (paleomagnetism).

D1.8 Typical values of susceptibility and remnant magnetization  Magnetic susceptibility susceptibility (k) in SI units

Sedimentary rocks Metamorphic rocks Granites Basalt/gabbro Ultramafic rocks

0.0005 0.0030 0.0050 0.0600 0.1200

I  /I  r  r  i = ratio of remnant to induced magnetization

0.01 0.1 1.0 10.0 -

Kearey Figure 7.5

 MJU March 2006

224 D2.1 The internal component of the Earth’s magnetic field

D2.1.1 Historical background

  Petrus Perigrinus (1269) described a floating compass and a vertically  pivoted compass (dip needle). He showed that the field of a magnet was dipolar  by studying a sample of lodestone. At that time the origin of the magnetic field was believed to be the Pole Star or a magnetic mountain located at the North Pole.

●

  In 1574 Gerhard Mercator realized that magnetic north and geographic  North did not coincide. This represented a discovery of the declination of the Earth’s magnetic field. He explained this problem by  placing multiple magnetic poles on his map.

●

 Measurements of the inclination of the Earth’s magnetic field were described in 15 76  by Robert Norman.

●

1

 William Gilbert developed these ideas and was the first to propose that the Earth itself was acting as a large magnet with a dipolar magnetic field (i.e. the magnetism did not originate in the stars or features on the surface at high latitude). He defined the magnetic poles as places where a dip needle would point vertically downwards.

●

  Maps of declination were refined, including the 1702 map drawn by Sir Edmund Halley, who proposed that two north magnetic poles and two south magnetic poles were needed to produce the observed field. Halley also observed that the declination was varying over time.

●

 Subsequent analyses required an ever increasing number of magnetic poles to describe the increasingly complex detail in the magnetic field that was determined by new observations. A more suitable mathematical description using spherical harmonics was developed by Gauss in 1839.

●

2

D2.1.2 Spatial variations

  The dipole accounts for about 80%  of the observed magnetic field. The axis of the dipole is inclined at  11.5 degrees  to the Earth’s rotation axis. The magnetic dipole 22 2 moment, P = 7.94 x 10  Am  in 2000.

●

 At any point the magnetic field is defined by the magnetic field elements

●

F = total total field field strength strength Z = vertical component of F H = horizontal component of F i = inclination (angle between F and surface) d = declination (angle between H and geog. north)

  The magnetic field is described by the  International Geomagnetic Reference Field  ( IGRF   IGRF ). ). The figure below shows the IGRF in 2000 and the database is updated regularly as the magnetic changes over time. Note the departures from a pure dipole field, especially the four regions of high F in high latitudes.

●

3

D2.1.3 Temporal variations in the internal magnetic field

  The compass was invented in China, and ancient texts record the variation of declination from at least AD720. It took a while for Europeans to catch on, and the first record of temporal changes was in 1635 when Gellibrand noted that declination varies with time. The declination in London was found to vary significantly over the period 1600-2000.

●

●

Secular variation in the Earth’s magnetic field occurs on many timescales including:

(a) Westward drift: features can be seen to move west over the last century.

4

(b) Investigations of historical records from early navigators and explorers has extended these records back to the 1600’s (Jackson et al., 2000) and are displayed as movies at http://geomag.usgs.gov/movies/ (c) Short term geomagnetic jerks occur on time scales of a decade

22

(d) Continuous reduction of dipole field  since 1600 ( P  = 9.4 x 10 22 2 1600, P = 7.94 x 10  Am  today)

2

Am in

(e) Complex sequence of magnetic field reversals  over the observed geological record. During a reversal the whole field switches north and south poles. Between reversals there is evidence that the magnetic dipole axis and the Earth’s rotation axis are approximately parallel. The sequence of reversals appears to be chaotic with no regular frequency. The present normal polarity (Brunhes chron) has lasted for 780,000 years.

5

7

8

Periods without a reversal for 10 -10  years are called superchrons. Cretaceous normal superchron 118-83 Ma Permo-Carboniferous (Kiaman) reverse superchron 312-262 Ma The reduction in the main field over the last 400 years could indicate we are approaching a reversal. How might reversals affect life on Earth?

D2.1.4 Origin of the internal magnetic field

  High temperatures inside the Earth (well above Curie temperature) exclude the  possibility of remnant magnetization generating the magnetic field. There is no large  bar magnet inside the Earth!

●

  The secular variation, and alignment of dipole with rotation axis, suggest that the magnetic field originates in the relatively rapid fluid motion in a part of the Earth with a

●

6

high electrical conductivity. This only leaves the outer core (composed of liquid iron) as the place where the magnetic field is generated.   A complex fluid motion is believed to act as a self sustaining dynamo. Computer simulations of the geodynamo  can partially explain the observed spatial and secular field variations, including reversals. These models include convection, Coriolis forces and magnetohydrodynamics. However many details remain unanswered.

●

Computer simulation of a geomagnetic reversal (Glatzmaier and Roberts, 1995)

D2.2 External component of the Earth’s magnetic field

The magnetic field measured at the surface of the Earth is due to sources inside  and outside  the Earth. The external component is generated in the atmosphere and magnetosphere.

The solar wind (a stream of H  and   and He ions) is deflected by the Earths internal magnetic field to create the magnetosphere.

●

7

The interactions between the solar wind and the Earth’s magnetic field are very complex. Temporal changes in the solar wind, due to sunspots, solar flares and coronal mass ejections can produce a change in the magnetic field at the surface of the Earth.

●

 From 50-1500 km above the Earth’s surface is the ionosphere, a region of plasma with high electrical conductivity. Changing magnetic fields from the magnetosphere can induce large electric currents in the ionosphere. Changes in these currents produce large changes in the magnetic field measured at the Earth’s surface.

●

 Large currents flow in specific locations including: -equatorial electrojet flows on magnetic equator on side facing sun -auroral electrojet flows at high magnetic latitude

●

  When the solar wind is in a steady state, the Earth’s magnetic field shows a daily variation that is due to the Earth turning within the current systems of the magnetosphere and ionosphere. The typical variation is called the solar quiet day

●

8

variation (Sq). The amplitude is typical 10-20 nT and varies with latitude. Clearly seen in time series above.

 A much smaller variation is seen every 25 days and is caused by the orbit of the moon.

●

  When the solar wind is active, the Earth’s magnetic field is said to be disturbed. Magnetic storms occur when the current systems change over a period of several days and the field at the Earth’s surface can change by 100’s of nanotesla. These changes are largest beneath major ionospheric current systems. A small substorm can be seen in the middle of the time series plotted above.

●

  Smaller magnetic field disturbances are classified as substorms and bays and have timescales of several hours.

●

  Solar activity is characterized by an 11 year cycle  and we have just passed the maximum. Maximum solar activity results in high levels of activity in the Earth’s external magnetic field and frequent magnetic storms and strong auroral displays.

●

References

G.A. Glatzmaier and P.H. Roberts, A three-dimensional self-consistent computer simulation of a geomagnetic field reversal, Nature, 377, 203-209 (1995). Jackson, A., Jonkers, A. R. T. & Walker, M. R., 2000. Four centuries of geomagnetic secular variation from historical records, records, Phil. Trans. R. Soc. London, A 358, 957-990.

9

D2.3 Comparison of the Earth’s Earth’s gravitational and magnetic fields

Gravitational field

Magnetic field

Overall field geometry

Approximate spherical symmetry g varies as 1/r 2

80% dipole B varies as 1/r 3

Direction

Down, by definition

Inclination varies from +90  to –90 ˚

˚

978,000 mgal at Equator 983,000 mgal at poles

25,000 nT at Equator 61,000 nT at high latitude

GRS formula simple and accounts for variation of g with latitude

IGRF is a complex series of spherical harmonics

Temporal variations with internal origin

Signal produced by plate motion and mantle convection????

Secular variation, jerks, westward drift and north-south field reversals Poles moving at ~ 15 km/yr

Temporal variations with external origin

Tidal signals (< 0.5 mgal)

Diurnal S q variation (50 nT) Magnetic storms (100-1000nT) 11 year sunspot cycle

Latitude variation in Edmonton

~ 1 mgal km -1

~3 nT km-1

Elevation variation in Edmonton

~ 0.3 mgal m-1

~ 0.03 nT m -1

Spatial variations

10

224D3: Magnetic anomalies produced by simple geological structures

Remember that objects can acquire both induced and remnant magnetization. Induced magnetization will disappear when the applied magnetic field is removed. Remnant magnetization is frozen into the material.

In the following examples, we will consider just induced magnetization. This will be in a direction parallel to the Earth’s magnetic field. Remnant magnetization can be in any direction.

3.1 Dike (Monopole)

the presence of the Earth’s magnetic field, the dike develops an induced magnetic moment

●  in

●  if

the structure extends to depth, then the lower monopole can be ignored since for a monopole, Br = μm/r 2

● compute

total field at surface by adding B r  and BE as vectors

3.2 Cylinder

cylinder has an induced magnetization with negative monopoles on upper surface and  positive monopoles on the lower surface. surface . One can show that net effect is a dipole at centre of cylinder 3 ● for a dipole, the magnetic field falls away as 1/r  ● note that the magnetic field anomaly is different at the magnetic north pole and equator ● The

3.3 Sphere ● As

for cylinder, the sphere develops an induced magnetic moment, which is equivalent to a dipole at the centre.

● what

happens in between? MATLAB demonstration

In map view

From http://www.gif.ubc.ca

3.4 Thin sheet or slab ●  an

infinite sheet develops negative poles (South poles) on the upper surface and positive  poles (North) poles on the lower surface. If the t he sheet is thin (vertical (vert ical direction), then these cancel to give no net magnetic field at the surface.

●  when

finite in horizontal distance, the poles near the edge do not cancel and a positivenegative anomaly is observed.

● thus

the magnetic anomaly is sensitive to the edges of structures

3.5 Basement topography ● superimpose

an infinite layer and a finite layer (D3.4)

 MJU March 2006 

224 D4 Instrumentation and data collection techniques for magnetic exploration D4.1 Instrumentation

D4.1.1 Flux gate magnetometer

 Details of operation are described in Keary and Brooks

●

 Measures the component of magnetic field along the axis of the ferrite cores

●

D4.1.2 Proton precession magnetometer

 cylinder contains a proton (hydrogen) rich liquid

●

 protons have a magnetic moment and act as small magnets

●

 they align with the Earths magnetic field, BE ● when a stronger magnetic field (BP) is applied, the protons are re-aligned ●

 when BP is switched off, the protons return to be being aligned with BE

●

 as this happens, they presses at a characteristic frequency, ω, that is proportional to BE

●

 this frequency is measured by the weak radio signal emitted by the precessing protons

●

 measures total magnetic field, so no alignment of the instruments is needed

●

D4.1.3 Alkali vapour magnetometer

 These magnetometers use changes in the frequency splitting of optical spectral lines of elements such as rubidium, caesium or potassium to measure the magnetic field. This Th is is governed by the Zeeman effect.

●

 Typically measure 10 times per second with precision of 0.01 nT or lower.

●

 Now the most widely used magnetometer in airborne surveys.

●

 Errors in survey limited by positioning errors rather than magnetic field measurements.

●

  Typical HRAM (High resolution aeromagnetic surveys) use these instruments. Elevation is 40-80 m and line spacing 200-500 m. This has allowed aeromagnetic exploration to be used in studies of the subtle magnetic anomalies generated by sedimentary rocks with small (but variable) magnetite contents. See examples in D6 from North Slope of Alaska, Canadian Cordillera and Weyburn, Saskatchewan. Example from Finland below.

●

D4.1.4 Magnetic gradiometers

 By placing one magnetometer above another, the vertical gradient of the magnetic field (dF/dz) can be measured. The high precision of modern instruments means the spacing only needs to be 1-2 m. This has several advantages

●

  By computing (dF/dz)/F, temporal variations in the external magnetic field can be effectively removed from the data. This is illustrated in the data example below (panel f), and this provides an alternative way to remove time variations (if the base station fails!)

●

  dF/dz can be more sensitive to dipole anomalies in the Earth than long wavelength features in the regional magnetic field.

●

D4.2 Field techniques for magnetic exploration D4.2.1 Ground-based magnetic surveys

 typically use a proton precession magnetometer ( ~ 1nT sensitivity) or an alkali vapor magnetometer (~ 0.01 nT sensitivity)

●

 operator must be liberated of all metal objects (especially their keys and money!)

●

 survey along a profile or grid

●

 locate profiles with surveying, or with an instrument with integrated GPS (above)

●

Spatial corrections

 It is possible to correct magnetic data for latitude and elevation but these corrections are generally much smaller than magnetic field anomalies. They are only important if the survey has a large spatial extent.

●

 Topography is important, such as when whe n measurements are made in a narrow valley va lley or canyon.

●

Time corrections

  A simple magnetic survey could use a single magnetometer that is carried along a  profile. Since the magnetic magne tic field of Earth changes with time, it is unclear if changes in magnetic field are due to (1) temporal changes or (2) subsurface structure.

●

To separate effects (1) and (2) two magnetometers must be used in fieldwork. An example of this procedure is illustrated below. One magnetometer is placed at a fixed location (the base station) and the other is carried along the profiles (b). Data recorded here are shown in red and measurements at the moving station are shown in blue (a).

●

  Variations at the base station instrument are subtracted from those on the moving instrument to determine variations along the profiles (c)

●

 Note that the three crossings of the target are virtually identical when the time variations are removed (d,f)

●

 Unlike gravity exploration, instrument drift is not a serious problem

●

D4.2.2 Aeromagnetic surveys

  Profiles Profiles flown with line line spacing 100 m to several kilometers. Elevation 1001000 m.

●

 Alkali vapour magnetometer. Accuracy 0.01 nT or less

●

 Continuous recording base station used to remove temporal variations.

●

  Cross lines used as internal check that temporal variations have been consistently removed.

●

  Need a non-magnetic aircraft or put the magnetic sensor in a towed bird. Compensation used to remove residual magnetic effects of aircraft.

●

 Flight path can be located with radar, altimeter or by taking aerial photos. GPS now makes this much more accurate and simpler.

●

Magnetic field data is usually corrected to a common elevation with upward or downward continuation

●

 Draped surveys flown at constant elevation above the surface (need helicopter)

●

 With total magnetic field measurements the motion of the aircraft or bird is not serious.

●

  Gives rapid coverage and works well in areas where access is difficult on ground (both in remote wilderness areas or urban environment)

●

  The magnetometer is further from target than in ground surveys, so the sensitivity is slightly lower. However, with alkali vapour magnetometers small ore deposits and features in sedimentary rocks can be mapped..

●

  High resolution aeromagnetic (HRAM) surveys are now the standard field technique and they use a dense line spacing  and low flight elevation  to give detailed information about subtle, sub-nanotesla magnetic field anomalies. See examples from Geological Survey of Finland at

●

http://www.gsf.fi/aerogeo/eng0.htm

D4.2.3 Marine magnetic surveys

  The proton precession magnetometer needs to be as far from a metal ship as  possible. This is usually in a plastic fish that is towed behind the ship.

●

 Marine surveying is slower and more expensive than airborne surveying, so not widely used in commercial offshore exploration.

●

  However, magnetometers are routinely towed behind survey ships and have revealed many details of the structure of ocean basins (e.g.  magnetic stripes  produced by seafloor spreading)

●

 MJU April 2006

Geophysics 224

D5 Interpretation of magnetic anomalies

5.1 Half-width techniques

  In section D4.1 we considered the magnetic anomaly due to a vertical dyke at a depth d  below the Earth’s surface (essentially a monopole). The half width was given by

●

x½ = 0.766 d  Other approximate rules can be derived but they must be used with considerable o f magnetic anomalies). caution (see Kearey p. 107 Interpretation of 5.2 Reduction to the pole

 Magnetic anomalies generally have a simpler shape when i= 90º than i = 0º

●

 Thus, interpretation would be easier if we could somehow change the inclination of the Earth’s magnetic field. This is difficult to do in practice.

●

  Reduction to pole uses mathematical filtering methodology. Kearey Figure 7.23 shows an example of synthetic magnetic data when i=60  and d = 20

●

˚

˚

  The measured data is shown in (b) with both positive and negative anomalies. Reduction to the pole (c) shows a simpler set of anomalies that reveal the true shape of the magnetized body.

●

5.3 Upward continuation  mathematical method for predicting magnetic field at a different elevation from that at which the data were collected.

●

 can allow data to be filtered, since at higher elevation the short wavelength features will be weaker than at ground level.

●

 synthetic example in Kearey figure 7.23(d) that was shown above.

●

5.4 Aeromagnetic maps

 Widely available and often produced by government agencies.

●

 Upward and downward continuation is needed to merge data collected at different elevations.

●

  Aeromagnetic data is often used for finding the lateral extent of structures or for locating anomalies that are then investigated in more detail during follow up on the ground.

●

 Anomalies are smaller when measured at flight elevation than at ground level. This can be both good and bad.

●

5.5 Two-dimensional and three-dimensional computer modelling

  Represent the Earth as series of complex polygons, each with a different susceptibility (k ) ●

 Can usually model both induced and remnant magnetization.

●

 both 2-D and 3-D computer codes are widely used

●

 Generally the user alters the model until the predicted and measured field data agree

●

 Automated inversion algorithms can be used

●

 Example of Winglink package Winglink  package with data from Alaska in D6.7

●

Geophysics 224 D6: Common applications of magnetic exploration methods D6.1 Detection of voids

 The air in a cave has k =0. If the host rock has a non-zero magnetic susceptibility, then a magnetic anomaly will be observed at the surface. What will be the sign and shape of the anomaly over a cylindrical lava tube in a basalt lava flow?

●

●

Ground-based magnetic surveys have been used in Hawaii to locate lava tubes prior to developing agriculture or construction.

●

Tunnels can be detected through their negative susceptibility contrast, as described above (case study at Teotihuacan in Mexico described by Arzate et al, 1990)

D6.2 Archaeology

●

●

Disturbing the soil can produce a small change in the magnetic susceptibility. This can permit the detection of foundations, graves or ditches (Clark, 1986)

Buried metal objects can be detected through induced or remnant magnetization

D6.3 Environmental geophysics ●

locating buried metal objects in landfills (55 gallon oil drums, pipes etc)

●

locating disused well casings. Why does the anomaly change with aircraft elevation?

  Mapping landfill boundaries (Roberts et al, 1990).

●

  Note that upward continuation  is a mathematical technique that computes the magnetic field at a higher elevation than that at which the data were collected.

●

  Why is the upward continued data smoother than the ground level data?

●

D6.4 Mineral exploration

 ore bodies can have a high magnetic susceptibility and may exhibit both induced and remnant magnetization. These anomalies can easily be detected at aircraft elevations and thus aeromagnetic exploration is a good reconnaissance tool for minerals.

●

 Example from Quebec (Kearey figure 7-25). The ore body is also electrically conductive, and is well defined by two geophysical surveys (electrical and magnetic) ●

  Remember that haemitite is antiferromagnetic and produces a negligible magnetic anomaly.

●

  Combined gravity and magnetic surveys sometimes help define if a magnetic anomaly represents a significant ore deposit. (Kearey figure 7-28).

●

  What is the physical basis of this technique?

●

 Examples from Atlin, British Columbia

Ultramafic body in Surprise Lake Valley P = peridotite (k = 88.5 x 10-3 SI) M = metabasalt k = 0.41 x 10-3 SI)

Mapping plutons and their associated mineralization. Previous geological mapping was revised. Magnetic field data courtesy of Carmel Lowe,  Natural Resources Canada.

P

D6.5 Diamond exploration Locating kimberlite pipes

Kimberlites are volcanic rocks that originate at depths of 100-200 km in the asthenosphere and move rapidly to the surface. If they originate in the diamond stability field, they can bring diamonds to the surface. Aeromagnetic data are widely used to locate kimberlite  pipes.

Kimberlite usually contains more magnetite than the host rock in the Slave Province and  produces a positive magnetic anomaly. However a negative magnetic anomaly occurs if: (a) the host rock is more magnetic (contains more magnetite) or (b) the pipe has a (reversed) remnant magnetization. If these effects cancel, the pipe can produce a weak , or near zero, anomaly. Extra information can be gained from airborne electromagnetic (EM) surveys (more details in Geophysics 424 next year). This is essentially a way of measuring shallow electrical resistivity from a moving aircraft or helicopter. Two factors make this a viable exploration technique: (1) Kimberlite has a lower electrical resistivity than the host rock (2) Kimberlites weather and produce clay, with a low resistivity. If eroded, a lake forms, also producing a lower resistivity. “Because kimberlite pipes exhibit variable anomalies on both electromagnetic and magnetic data, the best approach to mapping them is to simultaneously collect EM and magnetic data from a low flying platform. The DIGHEM system collects both data sets from a sensor at 30m altitude, sampling about every 3m.” 3m.” http://www.fugroairborne.com.au

(a) Point Lake Kimberlite, NWT

First kimberlite pipe discovered in NWT  by airborne magnetics and EM. Note the negative magnetic anomaly (top). There is also a strong response in the EM data (GEOTEM channel 7, off-time)

(b) Ekati Mine, Lac de Gras, NWT

The Ekati Diamond Mine is exploiting five economic kimberlite pipes in the Lac de Gras region of the NWT. The pipes are named Panda, Koala, Misery, Fox and Leslie. BHP Billiton mining operations at the Koala Pipe are shown on the right. Magnetic data: Left panel shows that the Grizzly pipe has a negative magnetic anomaly and Panda has a small positive anomaly. Koala and Fox pipes show weak anomalies. Airborne EM data: The apparent resistivity map (centre), calculated from the 7200Hz coplanar data of the DIGHEM survey, clearly shows the economic pipes in this data  block. The Fox pipe (south west corner) has the most distinct anomaly, and coincides almost exactly to the overlying lake. The Koala and Panda pipes give clear anomalies, and are also underneath lakes.

http://www.fugroairborne.com.au http://www.mining-technology.com/projects/ekati http://ekati.bhpbilliton.com

(c) Fort à la Corne kimberlites, Saskatchewan

Located beneath 100 m of sedimentary rocks and glacial overburden with no surface expression. Magnetic data are shown on the left and coincident GEOTEM data shown on right. www.fugroairborne.com

(d) James Bay Lowlands From Hogg and Munro (2000)

80 m elevation

30 m elevation

D6.6 Regional crustal structure Alberta Basement: The crystalline  basement rocks in Alberta date from the Archean and Proterozoic. However, they are covered by the sedimentary rocks of the Western Canada Sedimentary Basin and cannot  be studied directly. The basement rocks have been mapped through  potential field data (magnetic and gravity) and analysis of rocks recovered form the bottom of oil wells (Pilkington et al, 2000).

Even if the origin of the magnetization is not resolved, the character of the aeromagnetic anomalies can be used to determine the extent of a geological province. The direction of the magnetic anomalies can also reveal the geological strike of these rocks. Tibetan Plateau: Several geophysical techniques have suggested that unusually high crustal temperatures exist beneath the Tibetan Plateau. How will this alter the susceptibility of the crustal rocks? Good coverage with magnetic data in Tibet is hard to obtain on the ground (no roads) and aeromagnetic data coverage is not widely available.

In the 1990’s a low orbit satellite (MAGSAT) was used to map the Earth’s magnetic field. (Why in low orbit?) Analysis of these data by Alsdorf and Nelson (1999) reveal a  pronounced magnetic low over Tibet? Can this magnetic low be explained on the basis of high crustal temperatures and partial melting?

D6.7 Dike location Southern Alberta

These data were collected during the University of Alberta Geophysics field school in Southern Alberta. Shading shows magnetic field anomaly in nanteslas. Contour lines show surface elevations. The magnetic anomaly is due to a shallow dike.

Katmai Volcano, Alaska

The goal of the magnetic survey was to determine if a dike connects Mt. Katmai and  Novarupta. Both erupted in the 1912 1912 eruption, but the plumbing plumbing of this volcanic system is still unknown.

 Magnetic anomaly map:  map:  Note the eastwest magnetic high (approximately 1000 nT) that is coincident with the hypothesized dike trend. These data are the results of fieldwork by Graham Hill and John Eichelberger, Alaska Volcano Observatory, University of Alaska, Fairbanks.

 Model for western profile:  2-D model with variable magnetic susceptibility that fits data for one of the profiles.

Sedimentary basement  Pyroclastic flow  Rhyolite dike

k=0.004 k=0.002 k=0.01

D6.8 Hydrocarbon exploration with magnetic data

 A summary of the state-of-the-art in oil and gas exploration can be found in Gibson and Millegan, 1998. While oil and gas are not magnetic, useful information can be obtained from magnetic exploration, since it can define the local geology.

●

 Magnetic data can be used to define the depth of magnetic basement. The magnetic  basement is usually the crystalline basement and knowledge of the depth is useful in hydrocarbon exploration since it gives information about the overlying sedimentary rocks (depth, location of faults etc).

●

 Example: East China Sea, Okuma et al, in Gibson and Millegan, 1998, pages 59-62

  The Bowser Basin in British Columbia is currently being evaluated for hydrocarbon  potential. Note that spatial extent is defined from the magnetic data. Magnetic anomalies have longer wavelengths within the basin (why?). Magnetic field data courtesy of Carmel Lowe, Natural Resources Canada. ●

 Salt is diamagnetic and thus produces an anomaly of negative sign compared to a paramagnetic rock. This allows the geometry of salt diapirs to be defined from very accurate magnetic field data.

●

 Example: Gulf of Mexico Continental slope, Corine Prieto, in Gibson and Millegan, 1998, pages 14-16

 As oil seeps to the surface from a trap, it can alter the rocks through which it flows. This can change the near-surface magnetic susceptibility. Attempts to locate these altered regions have been made with magnetics, and other airborne geophysical data.

●

 Example: ALTREX method, J. Rowe et al, Gibson and Millegan, 1998, pages 124-129

 Note discussion of this effect by Machel and Burton (1991)   Generally sedimentary rocks have low magnetic susceptibilities and do not exhibit a strong induced or remnant magnetization. Sedimentary rocks can develop a weak remnant magnetization during deposition.

●

In detrital remnant magnetization  (DRM), magnetic mineral grains are oriented by the Earth’s magnetic field as they are deposited. With very sensitive magnetometers (alkali vapour) and accurate navigation in a high resolution aeromagnetic (HRAM) survey, the remnant and induced magnetization can be detected and interpreted.  Example: ANWR, North Slope of Alaska, Phillips et al, in Gibson and Millegan, 1998,  pages 130-134

 Example : Northern Canadian cordillera, shown in Nabighian et al., (2005)

In this case, the magnetic data can be used to map faults, in between the widely spaced seismic reflection lines.  Example: Weyburn carbon dioxide sequestration project, shown in Nabighian et al., (2005). The prescence a sedimentary layer with more magnetite than the background allows faults to be mapped in HRAM data. The seismic data are essential to ground truth the magnetic field data in this case.

 Example : Mushayandebvu and Davies, 2006. Magnetic gradients gradients in WCSB

D6.9 Seafloor magnetic anomalies

As new oceanic crust is formed at mid-ocean ridges, it is magnetized in the direction of the Earth’s magnetic field. For basalt, this remnant magnetization is usually stronger than the induced magnetization.

●

As the plates move away from the ridge, they carry a record of the magnetic field  polarity at the time the crust was formed.

●

These so-called magnetic stripes were a key piece of evidence in the development of  plate tectonic theory. By combining the history of the magnetic field reversals with radioactive dating, the age of oceanic crust can be determined from the magnetic field anomalies.

●

References

Alsdorf, D. and K.D. Nelson, Geology, 27, 943-946, 1999. Arzate, J.A., L. Flores, R. Chavez, L. Barba and L. Manzanilla, Magnetic prospecting for tunnels and caves in Teotihuacan, Mexico, in Geotechnical and Environmental Geophysics, Volume 3, SEG Investigations in Geophysics, No. 5, p. 1-30, 1990. Clark, A.J., Archaeological geophysics in Britain, Geophysics, 51, 1404-1413, 1986. Gibson, R.I., and P.S. Millegan, Geologic applications of gravity and magnetics: case histories, Society of Exploration Geophysics, 1998. Hogg, H., Munro, S., The aeromagnetic discovery of kimberlites and sulphides at depths up to 200m , Expanded Abstract SEG Annual Meeting, Calgary, 2000. Machel, H.G., E. Burton, Chemical and microbial processes causing anomalous magnetization in environments affected by hydrocarbon seepage, Geophysics, 56, 598-605, 1991. Macnae, J. C., Kimberlites and exploration geophysics, Geophysics, 44, 1395-1416, 1979. Mushayandebvu, M.F., and J. Davies, Magnetic gradients in sedimentary basins: Examples from the Western Canada Sedimentary Basin, The Leading Edge, January, 69-73, 2006. M. N. Nabighian, V. J. S. Grauch, R. O. Hansen, T. R. LaFehr, Y. Li, J. W. Peirce, J. D. Phillips, and M. E. Ruder, The historical development of the magnetic method in exploration, Geophysics, 70, 33ND-61ND, 2005. Pilkington, M, W.F. Miles, G.M. Ross and W.R. Roest, Potential field signature of buried PreCambrian basement in the Western Canada Sedimentary Basin, Canadian Journal of Earth Sciences, 37, 1453-1471, 2000. Roberts, R.L., W.J. Hinze and D.I. Leap, Data enhancement procedures on magnetic data from landfill investigations, in Geotechnical and environmental geophysics, Volume 2, Environmental and groundwater, SEG Investigations in Geophysics, No. 5, 261-266, 1990. Wilson, M. G. C., Diamonds Through the Decades: A Review of South African Production, Geotimes, July, 14-18, 1997.

 MJU March 2006

Geophysics 224 D: Magnetic exploration methods - summary

D1: Basic of geomagnetism 1.1

Definitions of Susceptibility (k) and permeability ( μ) Magnetic flux density ( B) and magnetic field strength ( H) Magnetic moment ( M) and intensity of magnetization ( I) 1.2 Monopoles 1.3 Magnetic field lines 1.4 Magnetic dipoles 1.5 Diamagnetism and paramagnetism 1.6 Ferromagnetism and domains 1.7 Curie temperature 1.8 Typical values of susceptibility and remnant magnetization

D2: The Earth’s magnetic field 2.1 Internal component of the Earth’s magnetic field 2.1.1 Historical background 2.1.2 Spatial variation 2.1.3 Temporal variation in the internal magnetic field 2.1.4 Origin of the internal magnetic field 2.2 External component of the Earth’s magnetic field 2.3 Comparison of the Earth’s gravitational and magnetic fields

D3: Magnetic fields of simple geological structures 3.1 Monopole 3.2 Dipole 3.3 Sphere and cylinder 3.4 Thin sheet 3.5 Basement topography

D4: Instrumentation and data collection techniques t echniques 4.1 Instrumentation 4.1.1 Flux gate magnetometer 4.1.2 Proton precession magnetometer 4.1.3 Alkali vapour magnetometer 4.1.4 Magnetic gradiometer 4.2 Field techniques 4.2.1 Ground magnetic surveys 4.2.2. Aeromagnetic surveys – HRAM surveys 4.2.3 Marine magnetic surveys

D5: Data interpretation 5.1 Half-maximum techniques 5.2 Reduction to pole 5.3 Upward continuation 5.4 Aeromagnetic maps 5.5 2-D and 3-D computer modelling

D6: Applications 6.1 Detecting voids 6.2 Archaeology 6.3 Environmental geophysics 6.4 Mineral exploration 6.5 Kimberlite exploration 6.6 Regional crustal structure 6.7 Dike location 6.8 Hydrocarbon exploration with magnetic data 6.9 Seafloor magnetic anomalies  MJU March 2006 

Geophysics 224 Winter 2006 Assignment 1 – Gravity methods 1. The LRT tunnel has a diameter of 3 m and is buried in sedimentary rocks with a density of 2600 kg m -3. The centre is at a depth of 4 m below the surface surface in downtown Edmonton. A gravity survey was run on a surface profile that crossed a long, straight section of the tunnel at right angles. (a) Sketch the variation in gz that would be observed. Calculate the minimum value of g z in mgal. Also calculate the half-width of the curve. (b) What is the minimum horizontal spacing  of gravity measurements that would be needed to detect the tunnel in (a)? Explain you answer carefully. (c) The gravity data has standard errors of 0.01 mgal. What is the maximum depth  at which the tunnel will be detectable on the basis of surface gravity measurements? Explain your answer.

2. Gravity data was collected at three stations (A,B and C) in Southern Alberta. These stations are on an East-West profile with latitude of 50  N. Station A is on the prairies, while B and C are further west in the Rocky Mountains. ˚

The absolute value of gravity at a nearby CSGN (Canadian Standardized Gravity  Network) station is 981,054.66 mgal. mgal. On the first day of the survey, the gravity crew make the measurements listed in the table below.

Location

Time

CGSN station Station A CGSN station

9:00 13:00 19:00

Relative gravimeter reading (mgal) 83.20 92.70 82.20

(a) Compute the difference in g  between the CGSN station and station A. This will require that you correct for the drift of the gravimeter. Plotting a graph of gravimeter reading as a function of time will help you do this. You may assume that the drift is linear with time. (b) Use your answer to compute the absolute value of gravity at station A. (c) What two factors cause the reading of a gravimeter to change when it is kept in a fixed location?

3. Gravity measurements were also made at stations B and C. The table below shows these measurements. Compute the Bouguer anomaly using the following steps. Please write your answers in the Table below. (a) Use the GRS67 formula from your notes to compute the predicted value of gravity (gnorm) at each station. Compute the gravity anomaly (g anom) by subtracting gnorm from gabs (b) Compute the Free air correction (CFA) for each station. (c) Apply the Free air correction to ganom to obtain the Free Air Anomaly at each station (d) Compute the Bouguer correction (C B) at each station, assuming a density of 2670 kg m -3. (e) Apply the Bouguer correction to obtain the Bouguer anomaly ( ΔgB) at each station. (f) What does the change in Bouguer anomaly from A to C tell us about the structure of the Rocky Mountains? Quantity Measured gravity gabs (mgals)  Normal gravity (GRS67) gnorm (θ = 50 ) Gravity anomaly ganom = gabs - gnorm Elevation above reference level (m) Free air correction CFA (mgals) Free air anomaly ΔgFA Bouguer correction CB (mgals) Bouguer anomaly ΔgB

Station A 981064.56

Station B 981002.93

Station C 980921.98

0

280

600

˚

You may assume that gravitational attraction of an infinite slab with thickness density  ρ is given by gz = 2πG ρΔ z and that the change in gravity due to a change in elevation of Δg = 0.3 Δ z (mgals)

 metres Δ z  metres

st

This assignment will be due in class on Tuesday January 31   2006 Please contact Ersan or myself if you have questions.

Δ z ,

is given by

Geophysics 224 Assignment 2 – Seismic Seismic reflection

1. The figure shows a shot gather that was recorded in a location where the subsurface consists of a series of horizontal layers. The travel times are listed on the attached table. You may answer the following by using graph paper or Excel. A copy of the travel times will be distributed by e-mail.

(a) For the first reflection plot a graph of normal moveout (NMO) as a function of 2  x ). Remember that the equation for normal moveout is an offset squared ( x approximation that is valid for  x << z1. Up to what value of  x is the approximation valid? (b) Compute values for the thickness ( z1) and velocity ( v1) of the first layer. 2

2

(c) To confirm your answer, plot t    versus  x   and compute the thickness ( z1) and velocity (v1) of the first layer. (d) For reflections 2,3 and 4, list the travel times at zero offset (t 2 , t 3 and t 4 ). 2

(e) Plot normal moveout as a function of  x   for reflections 2,3 and 4. Using the normal moveout equation

Δt  = t − t n =

 x

2

2 2t nV rms ,n

compute the r.m.s. velocities ( V rms, rms, 2, V rms, rms, 3 and  V rms, rms, 4). Be careful to only use travel times for which the equation is valid.

(f) Use the Dix equation to compute v2 , v3 and v4 1

⎡V rms,n t n − V rms,n −1t n −1 ⎤ 2 vn = ⎢ ⎥ t  t  − n n −1 ⎣⎢ ⎦⎥ 2

2

(g) From the travel times at normal incidence (t 2 , t 3 and t 4) compute the depths of the interfaces below the surface. th

(h) What does the shot gather tell us about the properties of the 5  layer? ______________________________________ x( m)

r ef 1( ms) r ef 2( ms) r ef 3( ms) r ef 4( ms)

______________________________________ - 200. 200. 000 000 235. 235. 36 477 477.. 26 723. 723. 17 966 966.. 07 - 190. 190. 000 000 233. 233. 08 476 476.. 55 722. 722. 84 965 965.. 89 - 180. 180. 000 000 230. 230. 87 475 475.. 86 722. 722. 53 965 965.. 72 - 170. 170. 000 000 228. 228. 76 475 475.. 21 722. 722. 24 965 965.. 56 - 160. 160. 000 000 226. 226. 75 474 474.. 59 721. 721. 96 965 965.. 41 - 150. 150. 000 000 224. 224. 85 474 474.. 01 721. 721. 69 965 965.. 27 - 140. 140. 000 000 223. 223. 06 473 473.. 47 721. 721. 45 965 965.. 13 - 130. 130. 000 000 221. 221. 37 472 472.. 96 721. 721. 22 965 965.. 01 - 120. 120. 000 000 219. 219. 80 472 472.. 50 721. 721. 01 964 964.. 89 - 110. 110. 000 000 218. 218. 35 472 472.. 06 720. 720. 81 964 964.. 79 - 100. 100. 000 000 217. 217. 01 471 471.. 67 720. 720. 63 964 964.. 69 - 90. 90. 000 000 215 215.. 80 471. 471. 31 720. 720. 47 964. 964. 60 - 80. 80. 000 000 214 214.. 70 470. 470. 99 720. 720. 33 964. 964. 53 - 70. 70. 000 000 213 213.. 73 470. 470. 71 720. 720. 20 964. 964. 46 - 60. 60. 000 000 212 212.. 89 470. 470. 47 720. 720. 09 964. 964. 40 - 50. 50. 000 000 212 212.. 17 470. 470. 26 720. 720. 00 964. 964. 35 - 40. 40. 000 000 211 211.. 58 470. 470. 09 719. 719. 92 964. 964. 30 - 30. 30. 000 000 211 211.. 12 469. 469. 96 719. 719. 86 964. 964. 27 - 20. 20. 000 000 210 210.. 79 469. 469. 86 719. 719. 82 964. 964. 25 - 10. 10. 000 000 210 210.. 59 469. 469. 80 719. 719. 79 964. 964. 23 0. 000 000 210. 210. 53 469 469.. 79 719 719.. 79 964 964.. 23 10. 10. 000 000 210. 210. 59 469 469.. 80 719 719.. 79 964 964.. 23 20. 20. 000 000 210. 210. 79 469 469.. 86 719 719.. 82 964 964.. 25 30. 30. 000 000 211. 211. 12 469 469.. 95 719 719.. 86 964 964.. 27 40. 40. 000 000 211. 211. 57 470 470.. 09 719 719.. 92 964 964.. 30 50. 50. 000 000 212. 212. 16 470 470.. 26 720 720.. 00 964 964.. 35 60. 60. 000 000 212. 212. 88 470 470.. 46 720 720.. 09 964 964.. 40 70. 70. 000 000 213. 213. 72 470 470.. 71 720 720.. 20 964 964.. 46 80. 80. 000 000 214. 214. 69 470 470.. 99 720 720.. 33 964 964.. 52 90. 90. 000 000 215. 215. 78 471 471.. 31 720 720.. 47 964 964.. 60 100. 100. 000 000 217. 217. 00 471. 471. 67 720. 720. 63 964. 964. 69 110. 110. 000 000 218. 218. 33 472. 472. 06 720. 720. 81 964. 964. 79 120. 120. 000 000 219. 219. 79 472. 472. 49 721. 721. 01 964. 964. 89 130. 130. 000 000 221. 221. 36 472. 472. 96 721. 721. 22 965. 965. 01 140. 140. 000 000 223. 223. 04 473. 473. 47 721. 721. 45 965. 965. 13 150. 150. 000 000 224. 224. 83 474. 474. 01 721. 721. 69 965. 965. 27 160. 160. 000 000 226. 226. 73 474. 474. 59 721. 721. 95 965. 965. 41 170. 170. 000 000 228. 228. 74 475. 475. 20 722. 722. 23 965. 965. 56 180. 180. 000 000 230. 230. 85 475. 475. 85 722. 722. 53 965. 965. 72 190. 190. 000 000 233. 233. 06 476. 476. 54 722. 722. 84 965. 965. 89 200. 200. 000 000 235. 235. 36 477. 477. 26 723. 723. 17 966. 966. 07 _______________________________________________________

2. Read chapters 3 and 4 in the textbook  th

This assignment will be due at class on Thursday March 16  2006. Please contact Ersan or myself if you have questions or simply need some inspiration to get started.

Geophysics 224 Assignment 3 – Seismic reflection reflection and refraction Question 1

The figure above shows an unreversed seismic refraction profile. (a) Compute the velocity at which the direct wave travels. Compute maximum and minimum velocities consistent with the measured data (i.e. estimate an uncertainty in v 1). (b) Compute the velocity (v 2) for the refracted arrival. Again, compute maximum and minimum velocities consistent with data (i.e. estimate an uncertainty in v 2). (c) What is the depth of the layer with velocity v 2? What are the maximum and minimum values that are possible? (d) A third layer is present in the study area with v 3 = 3000 m/s, but is not observed on the seismogram above. What is the minimum depth  at which it could be located? (e) What assumptions have you made answering (a)-(d)? Question 2

Clearly explain what is meant by : (a) interval velocity (b) r.m.s. velocity (c) stacking velocity.

Question 3

A 48-channel streamer has hydrophones spaced 10 m apart. As the survey ship moves forward, a shot is fired every 40 m. (a) What is the fold of CMP coverage? (b) To obtain 24-fold coverage, what must the new shot interval be? Illustrate each answer with a sketch

Question 4

The figure above shows two depth sections through the subsurface. Sketch the shape of the zero-offset time section that would be recorded by a profile that crosses the area. There is no need to compute exact travel times, but include the effects of diffractions etc in a qualitative manner. Velocities are in m/s and show relative velocities.

Question 5

Read chapter 5 from Kearey and Brooks, Seismic refraction methods

This assignment will be due in class Tuesday April 4th 2006.

MJU March 2006