Application of the Rayleigh-FFT technique to magnetotelluric modeling and correction

Physics of the Earth and Planetary Interiors, 53 (1989) 365—375 Elsevier Science Publishers B.V., Amsterdam — Printed in The Netherlands

365

Application of the Rayleigh—FFT technique to magnetotelluric modeling and correction G.R. Jiracek, R.P. Reddig

*

and R.K. Kojima

* *

Department of Geological Sciences, San Diego State University, San Diego, CA 92182 (USA.) (Received December 30, 1986; revision accepted August 31, 1987)

Jiracek, G.R., Reddig, R.P. and Kojima, R.K., 1989. Application of the Rayleigh—FFT technique to magnetotelluric modeling and correction. Phys. Earth Planet. Inter., 53: 365—375. A modification of the original 1896 Rayleigh scattering theory has been implemented using the fast Fourier transform (FFT) to calculate magnetotelluric (MT) responses for complex Earth structures including topography. The method has been tested for three-dimensional Earth models with arbitrary irregular layering and topography, although results are presented for two-dimensional (2-D) Earth models only. As layers can be pinched-out, isolated surface and sub-surface bodies can be included. An inherent error in the Rayleigh method limits the maximum surface and sub-surface slopes that can be accurately modeled. Valid 2-D results are obtained for surface slopes of 53° in the TE mode and 260 in the TM case; however, sub-surface interface slopes can exceed 60 ° for both modes. The Rayleigh—FFT approach enables simultaneous calculation of many surface points, typically 32 or 64 in 2-D. Model input requires only the digitized values of surface and sub-surface interfaces and the electrical properties of the various media. The point TM topographic distortions are slope dependent and are nearly height independent, e.g., maximum TM distortions for topography of 1 m and 1 km are nearly equal if the surface slopes are equal. However, since measured electric fields are averaged over finite length dipoles, the point-by-point small-scale effects are greatly smoothed in practice. TE topographic distortions are typical frequency-dependent inductive effects which are much smaller than TM effects as manifested in the apparent resistivity calculations. Use of the Rayleigh—FFT algorithm to remove the unwanted effects of topography and near-surface inhomogeneities is demonstrated by 2-D distortion tensor stripping of both theoretical and actual field MT results.

1. Introduction The magnetotelluric (MT) effects of surface topography and near-surface inhomogeneities, such as irregular sediment-filled basins or lava flows, can greatly exceed, and therefore mask, the MT response of sub-surface targets. Such effects result in a parallel displacement of the logarithmic apparent resistivity sounding curves at periods when the skin depth exceeds the size (vertical and lateral extent) of the distorting feature. This period *

**

Present address: Sohio Petroleum Company, Dallas, 1’X, U.S.A. Present address: Earth Sciences Associates, Palo ~to, CA, U.S.A.

0031-9201/89/$03.50

© 1989 Elsevier Science Publishers B.V.

range often begins at frequencies higher than those recorded in field data. In this case, the results at the shortest periods are already displaced, as is evident in the MT sounding in Fig. 1. The transverse electric (TE) and transverse magnetic (TM) values in Fig. 1 are shifted approximately parallel to each other from the beginning to 0.2-s period. Above 0.2 s, larger lateral Earth effects destroy the parallelism caused by the smaller-scale features. Distortion effects are termed ‘static’ in the period range where they produce a constant shift in logarithmic apparent sounding curves. This shift is caused by resistivity electric charge build-up or galvanic effects (Andrieux and Wightman, 1984; Reddig and Jiracek, 1984) at lateral resistivity

366

~iooo

]-l-rm

i

I

I

I

-

GRANDE

ioo

—

0

•

f

00 00

~

0 0

10

~

o

0

0 0

•.

-

J

U, —

_____

_____

_____

_____

-

OBSERVED DATA

w ~

~ 0.’ .LLIJII• 0.02 0.1

I

I I

iiii

I

11111

I

I I

0

iii

100

I

I

i

-

1000 2000

PERIOD (SEC) Fig. 1. TE and TM apparent resistivity sounding data at station 6-79 near the Rio Grande nft showing short- and long-period static shifts (after Jiracek et al., 1987).

changes, i.e., topography or inhomogeneities. The shift is frequency dependent if the normal component of the charging electric field is not uniform over the resistivity change. This is the case at frequencies high enough that the skin depth is small compared with the distorting feature. If the resistivity contrasts are two-dimensional (2-D), then only the true TM curves are shifted significantly. However, small-scale features are more likely to be three-dimensional (3-D); therefore, both apparent resistivity curves are shifted, usually by different amounts. This is so since both curves have contributions from electric fields norma! to boundaries in the 3-D case, i.e., there is no true TE—TM separation. A constant vertical shift of a single curve cannot correct the 3-D situation but both curves can be successfully shifted if the correct surface resistivity section is known. Both Andrieux and Wightman (1984) and Sternberg et a!. (1985) advocated auxiliary electromagnetic soundings to facilitate this correction. Jones (1988) suggested the shifting of all soundings within a basin to a modal or uniform surface resistivity. Another method is practiced in the U.S.S.R. where, in a region of nearly conformal but shifted sounding curves, all of the TE and TM curves are each logarithmically averaged to yield a single representative TE and TM curve for interpretation, e.g., Berdichevsky et al. (1980).

A completely different technique used in the U.S.S.R. to correct for surface inhomogeneities is to shift the longest period branches of the sounding curves into coincidence (Vanyan, 1981). This method (Rokityansky, 1982, pp. 201—209) assumes that, at depths greater than several hundred kilometers, the resistivity within a tectonic province does not vary laterally and that a long-period standard or normal apparent resistivity curve applies. Global geomagnetic depth sounding (GDS) daily-period data are sometimes used for tying the MT results. The long-penod method descnbed above addresses larger-scale static problems as produced by a large irregular sedimentary basin. Such a .

.

.

.

static shift is evident above 100-s penod in Fig. 1. Here, the TE and TM curves have reached another static shift caused by the galvanic effects at lateral boundaries some tens of kilometers away. For MT researchers interested in midcrustal or deeper targets, such long-period shifts represent distortion just as troublesome as the shorter-period static shifts. Clearly, it is desirable to correct for distortion effects of all scales. Curve shifting will not solve for these effects if statistical methods are employed when the area is inadequately sampled or if shifts are applied where distortions are not frequency independent (not static). Reddig (1984) showed that in the latter situation, false interpretations will result after leading the interpreter to believe that distortions have been removed. A new magnetotelluric field procedure and data-processing method called EMAP (Bostick, 1986) promises to counteract all undesirable distortion effects in 3-D and in 2-D TM data. These effects are removed by spatial low pass filtering via both continuous electric field (dipole) sensing and frequency-dependent, weighted averaging along a profile. Field examples of the EMAP procedure (Shoemaker et a!., 1986; Word et a!., 1986) are very encouraging but continuous profiling is laborious and it may be logistically untenable in some circumstances. Also, spatial filtering of the data may compromise vertical and horizontal resolution to some degree. Some success in simply numerically modeling the surficial and near-surface distortions along with the deeper targets has been achieved in 2-D

367

(Wannamaker, 1983; Jiracek et a!., 1987). Only simple block models have been studied in 3-D (Park et a!., 1983; Wannamaker et a!., 1984; Park, 1985). However, quasi 3-D thin sheet conductance calculations have proved useful (Kaikkonen et a!., 1985). A distortion tensor stripping technique as proposed by Larsen (1977) was used by Mozley (1982) in 3-D to correct approximately for the long-period topographic effects on MT data collected in mountainous terrain in the western United States. Tensor stripping of topographic effects is in the realm of possibility since topography is deterministic. Distortion tensor removal of near-surface inhomogeneities, however, is risky since the distorting causes must be completely defined or assumed. This paper describes an accurate, efficient numerical solution to many MT topographic problems in 2-D and 3-D. The method is based on the original Rayleigh (1896) hypothesis that all scattered fields from a rough surface can be expressed as superpositions of outgoing plane waves. The method has been extended to Earth models

with arbitrary numbers of irregular layers. Thus, distortion tensor stripping of the effects of known sub-surface inhomogeneities can be made with the Rayleigh modeling method.

2. Rayleigh scattering approximation Figure 2 illustrates Rayleigh scattering from a two-dimensionally rough surface, ~(x, y). This is a 3-D problem since the electromagnetic properties vary in 3-D. A plane wave source is assumed to impinge upon the rough surface at an arbitrary incident angle 9~ with the z-axis. The Rayleigh hypothesis leads to an expression (Jiracek, 1972, 1973) for the total (t) electric field, in the y direction, in medium 0 of Et (x ~ =

exp( 1 +

—

ik0x sin 0~+ ik0z cos 9~) ~ ±oc

~ j

~.

j

±

~ (kr, k~)

xexp(ik~x+ ik~y ik~0z)dk~dk~ —

(1)

INCIDENT PLANE WAVE FRONT

-t

HQ

FREE SPACE MULTIPLE SCATTERING

60

E0,/Lo,p0~OO fl-M

—S ~—S

E0,H0

m

\\

HALF SPACE

x

.—s —s E1, H1

Fig. 2. Schematic representation of Rayleigh electromagnetic scattering along a single profile of a two-dimensionally rough surface (after Jiracek, 1973).

368

Here, k~,k~and k~0are wavenumbers in the x-, y-, z-directions, and 2, (2) k~0= i(k~+ k~ k~)~’ 1/2 = w (~o ia 0/w )~ —

—

are designated Ao(k~,k~)and C0(k~,kr), respectively. The complete unknown A0, B0, and C0 coefficient sets act as reflection coefficients for the scattered x, y, and z electric field wavelets. These

—

As usual, w is the radial frequency; Ec~,a0, and Po are the free space permittivity, conductivity, and magnetic permeability. Since free space is considered to be a perfect insulator, the resistivity Po = 1/a0 = no, and k0 in (2) is a real quantity. Equations similar to (1), but without the initial sourceterm, are written for the x and z components of the electric field in medium 0 (Jiracek, 1972). The first term in (1) is that of the incident y-directed electric field propagating with direction cosines k~= k0 sin 0~and k~= k0 cos O~.The second term is the y-directed scattered field which, by the Rayleigh hypothesis, is a superposition of plane waves; the function B0(k~,k~)exp(ik~x+ ik~y 1k20z) (3) —

is an elemental plane wave where the sign of the radical k~0in (2) has been chosen to assure that the scattered field obeys the radiation condition as z + 00. That is, the imaginary part of k20 must be 0. The integration over all k~and k~in eqn. (1) (—no k~,k~ + no) results in two distinct types of elemental planeradiating wave contributions to the scattered field, namely and nonradiating terms. From expression (2), when —~

(k~+ k~)
(4)

quantities constitute the solution for the scattered electric fields above the rough surface. Scattered fields below the surface, in medium 1, are similarly written as superpositions of plane wavelets with x, y, and z electric field coefficients A1, B1, and C1( k~,kr). The radiation condition below the surface requires the z exponential term in the double integral (as in eqn. (1)) to be exp(+ik21z). Since the sub-surface is attenuating (i.e., a * 0), k1 is always complex as can be seen by substituting subscript 1 for 0 in (2). Now the scattered wavelets in medium 1 (and subsequent layers) will always have both a radiating and attenuating cornponent in z. These radiating wavelets are usually inhomogeneous waves since the planes of constant phase and constant amplitude typically do not coincide (Stratton, 1941, pp. 504—505). The basic Rayleigh assumption that the scattered fields are composed of plane wavelets that either propagate and/or are attenuated away from the rough surface went unchallenged for over half a century. It was then pointed out by Lippmann (1953) that above the surface, within the m) maximum (~(x, y)M)(Fig. and2),minimum surface excursions multiple (~(x,y) scattering should allow both upgoing and downgoing radiatrng waves (as well as exponentially growing and attenuating non-radiating waves). The Rayleigh

k~ 0is positive, real, and the wavelet (3) propagates away (upward) from the surface without loss. This is a radiating wave. When 2 (5) (k~+ kfl > ko k~ 0is negative, imaginary, and from expression (3) it is clear that there is only attenuation in the + z direction. There is propagation in the x and y directions but not in the z direction. These latter wavelets are non-radiating waves; they are also called evanescent, inhomogeneous, and surface waves. These non-radiating terms make an important contribution to the scattered electromagnetic fields on the surface. The coefficient sets for the scattered electric fields in the x and z directions above the surface

solution allows only upgoing radiating waves above the surface and downgoing ones below the surface. These are schematically represented by the solid wavelets that2; are not included Rayleigh rays in Fig. the dashed rays inin Fig.the2 represent hypothesis. The neglected terms in the Rayleigh method constitute an inherent error; however, in many instances this error is insignificant. As discussed below, the Rayleigh hypothesis is valid when the surface slopes are modest but sub-surface interface slopes can be quite large. 3. Rayleigh—FFT solution To evaluate the scattered electric fields, we must first solve for the six unknown coefficient

369

sets: A0, B0, and C0(k~,k~)above the surface and A1, B1, and C1 (kr, k~)below the surface. To do this we impose the boundary conditions on the continuity of tangential electric and magnetic fields. This provides four equations; the remaining two equations may be obtained by the continuity of the normal components of the magnetic induction and electric displacement or by setting the divergence of the total electric fields in each media (0 and 1) equal to zero. In either case, this yields a system of equations 4 j +~ j +~ ~A~g + f I I jl p 1g + 0g + 1g .

.

5+ C 6) C0g~ 1g~ exp(ik x + ik y) dk~dk

+

~

exp(—ikox sin 0~)

electrical properties the media. No grid mesh or fimte elements are ofneeded as required by dif-

j= 1,2,3,4,5,6 (6)

where ‘11k

1’ =

g

g

~

X~

k ~

x ‘

and ~

was first applied to electromagnetic scattering by Jiracek (1972, 1973) where he designated the approach the Rayleigh—FFT method. The Rayleigh—FFT approach enables the simultaneous calculation of electromagnetic scattering at many equally spaced surface points. This is typically 32 in both x and y directions for the 3-D problem or 32 or 64 in 2-D, i.e., a power of two as dictated by the standard FFT. The maximum number of discrete scattering wavenumbers (eqn. (1)) is always one less than the number of digitized surface points, e.g., 31 for both k and in 3-D. The model input for the Ray!eigh—FFT method (and requires the digitized valuesand of the surface any only interface) irregularities

ferencing schemes. Kojima (1985) extended the 2-D Rayleigh—FFT algorithm to an arbitrary number of layers. Since the layers can be pinchedout, isolated sub-surface bodies can be modeled. Rayleigh—FFT results for the 2-D problem have been extensively tested by comparison with other scattering theories, e.g., with integral equations

—

~X~ ~1

are known quantities dependent on the deterministic surface, ~(x, y). Detailed discussion of the above procedure may be found in Jiracek (1972, 1973) and Reddig (1984). A practical solution for the required coefficients is realized by first assuming a periodic rough surface. The infinite integration limits in (6) become infinite summation limits. These limits are subsequently truncated to finite values. The boundary conditions are then applied at equally spaced points along the surface (a point matching procedure) where the number of points in the x and y directions is greater than the truncated, finite number of the discrete k~and k~values. At this stage of the solution, a technique first used by Aki and Lamer (1970) in seismic scattering is used to solve the system of simultaneous, linear equations in the spatial frequency domain. The equations are better behaved in the wavenumber domain, and the final truncation to produce a square system can be thought of as data compression (Jiracek, 1972, 1973). The Aki—Larner procedure

and physical finite differenceoptics modeling (Jiracek, (Reddig, 1972, 1984), 1973), and with finite element results (Kojima, 1985; Wannamaker et a!., 1986). For perfectly conducting sinusoidal surfaces, the Rayleigh—FFT TE results are valid for maximum surface slopes of 310 (Jiracek, 1972, 1973); this is independent of the height of the surface topography except when the incident wavelength is much larger than the surface height. For very long wavelengths, such as those in magnetotellurics, the slope limitation has been defined by Reddig (1984) as a maximum of and 260 for the TE and TM modes, respectively. Irregular sub-surface interfaces with slopes exceeding 600 have been correctly modeled for both modes in the 2-D case (Kojima, 1985). Tests have not been completed for 3-D Rayleigh—FFT modeling but it is reasonable to expect that surface slope limitations of near 260 (similar to the 2-D TM mode) will apply. In all cases, the slope limitations are not as restrictive when the Earth mode! is more resistive. The most diagnostic tests of the validity of the Rayleigh—FFT results are derived from self-con—

—

530

370

sistency checks. These evaluate the degree to which the surface boundary conditions are matched and the amount of change in the results for increased numbers of scattering orders (truncations). These are expressed point by point or by overall rootmean-square values; the latter are typically less than a few per cent for valid results. We have found that in situations when the inherent Rayleigh error approaches an unacceptable level, the solution series often converges but then diverges as more scattering orders are used. Reddig and Jiracek (1984) have presented Rayleigh—FFT results to characterize the TE and TM topographic responses from 2-D corrugated (sinusoidal) surfaces. The maximum TM apparent resistivity distortion is almost entirely sensitive to the maximum surface slope. This is understood by realizing that the surface charges (galvanic effects) that produce the secondary fields responsible for the TM effects are greater when the assumed horizontal, incident electric field is more nearly perpendicular to the surface facets, i.e., the surface slopes are larger. Surface slopes with very low amplitudes (say 1 m) produce TM topographic distortion as large as surfaces with larger amplitudes (say 1 km) if the surface has the same maximum slope. However, the effects are spread out more for larger features and since, in practice, we measure electric fields with finite length dipoles, the effects of small features tend to be averaged out That is in reality we spatially low pass the data rather than make point measurements. The apparent resistivity distortion for the TE mode is much less than the TM case. The TE response pattern is more frequency dependent since it is an induction phenomenon.

4. Distortion tensor stripping Having established the applicability of the Rayleigh—FFT technique in MT modeling, we use it in a forward sense to remove the topographic and near-surface effects by distortion tensor stripping. We have extended the earlier work of Larsen (1977) and Moz!ey (1982) by allowing both the electric and magnetic fields to be linearly related to the undistorted fields by complex, frequency

dependent, three by three distortion [D] tensors. For example, the expression for this relationship between the x, y, and z components of the distorted (E’~)and undistorted (E”) electric fields is generally ED

DE

x

xx

E~ ED

=

z

DE

DE

EU

xy

D~ D~ D~

E~-’

DE

DE

EU

zx

ZY

DE

(7)

In practice, the undistorted Earth model is taken to be a flat, homogeneous half-space. Therefore, as z —* 0, the value of the undistorted vertical electric field (Ei’), measured in the Earth at the surface, vanishes since there is no normal component of current density. Also, a homogeneous half-space does not produce a vertical component of magnetic field (HZU) so this term is also equal to zero. The results presented in this paper are for 2-D geoelectric models where the y-axis is the invariant direction. Here, E~is the TM electric field and E~the TE case. Expression (7) (along with the magnetic equivalent) in the TM case now reduces to —

D

E

U

E~ ~ E~ D~E~’ HD DHHU ~ ~ and in the TE mode =

=

(8)

=

~

.v~ ~v

H~ D~HJ~’ HD DHHU =

(9)

=

Z

To implement eqns. (8) and (9) we calculate the distortion tensor elements by simply computing the complex ratios of appropriate distorted to undistorted fields. The distorted E and H fields are calculated by the Rayleigh— FFT method, which requires that we know the geometry and electrical properties of the distorting features. The undistorted fields are calculated assuming a homogeneous half-space as discussed above. We, of course, assume no a priori knowledge of the subsurface targets; therefore, distortion tensor stripping neglects electromagnetic coupling between the distorting features and the sub-surface targets.

371

5. Tensor stripping results

.2

We shall demonstrate the distortion tensor stripping method in conjunction with Ray!eigh—FFT modeling by considering two theoretical models and one actual field result. Figure 3 contains a 2-D model where the topography is cosinusoida! with 200-m relief and a wavelength of 3000 m. The bottom boundary of the 10 ~ m top layer has a graben-like depression from 1000 to 1400 m elevation. It is centered under the topographic valley. The periodicity mequirement in the Rayleigh—FFT approach results in the sub-surface interface actually being a horst—graben sequence. The lower half-space is a resistive basement of 100 m. Distortion tensor stripping has been applied by Reddig (1984) to synthetic MT data at. many sites along the surface, sites 1, 2,... (Fig. 3). Here, we will only present results at station 1 (x 0) on the topographic hill. Figure 4 contains plots of the real parts of the TE and TM distortion tensor elements versus period for the horizontal electric and magnetic fields (eqns. (8) and (9)). The value of the tensor elements is equal to 1 when there is no distortion, These plots contain the key information to understand the physics of the topographic distortion process. Both the electric and magnetic fields are distorted in the TE case but the amounts decrease —

—

~

=

I

I

I

I

I

I

I

I

I

I

I

I

£

MT SITES

~ O

-500

— -

~

o ~

1.0

TM

I— Z

~

—

0 I-

-

0 0.8

-

~

I

0,01

0.1

I

I 10

I 00

000

PERIOD (SEC) Fig. 4. Real parts of the TE and TM distortion tensor elements versus period at Site 1 (x = 0) in Fig. 3.

with increasing period (D~ and D~—s 1) as expected for an induction process. The TE distortion also approaches zero at much higher frequencies when the skin depth in the 10 ~2m layer is much less than the topographic relief. It is interesting to note that although both the TE electric and magnetic fields are distorted by 2-D topography, the response as manifested by the apparent resistivity distortion is small (Reddig and Jiracek, 1984). The fields are distorted in such a way (Fig. 4) that their ratio (Ef1/Hf), which is used to calculate the apparent resistivity, is little affected. The TM electric field distortion approaches zero (D~ 1) at short periods because of reduced skin depth. (This isslopes.) not strictly true for MT sites located on surface At longer periods, the skin depth is much greater than the surface topography and the surface charges along the topographic slopes increase and saturate. This galvanic effect reaches a near-constant value of distortion

10 A-M -

2—D MODEL VERT. EXAG. -1.33X

-

Ui

~

—.

.100

-

1Lx

~(X)0t~cos(21.rx/ x0

-

above period 4). This‘static’ asymptotic behavior is100-s the basis for(Fig. the constant shift —

-1000

1 -1500

-

0

100 fl~M\ I

I

I

//~o I

1000 2000 DISTANCE CM)

n-u I 3000

Fig. 3. Two-dimensional geoclectric model with cosinusoidal topography and horst—graben sub-surface feature (after Reddig, 1984).

as discussed above in the introduction. Although there is some magnetic field distortion in the TM case it is a very small fraction of the undistorted magnetic field and the resultant ratio D~ =

H~/H~(expression (8)) is essentially 1 at all frequencies (Fig. 4). Therefore, the 2-D TM topographic distortion is controlled solely by the elec-

372

90

-

~

80

-

0

70

-

E

60

-

(1)

50

—

~

40

-

—

I

I

100

I

SITE 1

-

TE

KEY

(0 & uD),/~

0-DISTORTED UD

/

UNDISTORTED

-

Figure 6 is a more complex theoretical, 2-D

I

/

-

TM

~IlS~:0RTIO;~ll

-

TM(D)

0.

~ ic 0.01

Figure 7 presents the appropriate model calcu-

-

DC)

3020 -

0

400-m-deep, 2042-rnThe basin which represent distorting features. target in this model two is a 1042—rn 7.5to 2.5-km crustaldepth. conductor with a bulge rising from

CORRECTED

Ui

Z

geoelectric model with a 600-m-high hill and a

I

I

I

0.1

1

10

100

1000

PERIOD (SEC)

Fig. 5. Distorted, undistorted, and distortion-corrected TE and TM apparent resistivity versus period curves at Site 1 in Fig. 3.

tric field distortion. The corresponding imaginary parts of the distortion tensor elements plotted in Fig. 4 are all <0.06 magnitude. They have a peaked character, like the imaginary electromagnetic response of isolated bodies (e.g., Ward, 1967, pp. 74—86), but their contributions to the distortion tensor elements appropriate to the model in Fig. 3 are insignificant. Figure 5 presents distorted (D), undistorted (UD), and distortion-corrected (DC) TE and TM synthetic sounding curves at site 1 in Fig. 3. The vertical scale is linear for clarity rather than the usual logarithmic scale. The TE distortion at site 1 is insignificant, hence, the TE (D) and (UD) curves are indistinguishable. However, the TM distortion results in a measured, anomalous decrease in apparent peaks. of The frequency TM (D) since curveresistivity D~ in Fig. is a 5on function is topographic clearlyofa period function in Fig. 4. Therefore, a ‘static’ shift would not properly correct the TM results at periods less than 100 s. Applying the Rayleigh—FFT calculated tensor corrections to the distorted TM (D) results produces the distortion-corrected TM (DC) curve in Fig. 5. This agrees with the TM undistorted (UD) results, i.e., a sounding measured as if site 1 was located on a flat surface of zero elevation,

to test distortion-correction stripping of the mode! in Fig. via lations 6. The the results Rayleigh—FFT are at one which period are(10 required s) for the TM mode. To calculate the distorted electric and magnetic fields (expression (8)) we first determine the response of a model with only the hill and basin and no deep crustal conductor. This result is designated ‘hill and basin’ in Fig. 7. It is combined with the response of a flat, homogeneous half-space of 100 ~2 m to correct the synthetic response of the complete model in Fig. 6 (the ‘hill, basin, and bulge’ curve in Fig. 7). These latter results are then tensor distortion-corrected to remove the effects of the hiJi and basin; the 17 .1

I

I

I I

I

I

HILL

-

I

I

BASIN

-

0 -

-i

—

—

100

~

n-u

-

—

2

—

z 0

~ -

> Ui .~l

ion-u

Ui

-s

BULGE

..6

-—

—7

—

2-D MODEL VERT. EXAG.

—

0

1I

2I

3I

4

5

6

7

8

9

10

DISTANCE (KM) Fig. 6. Two-dimensional geoelectric model with distorting hill and basin features above a crustal conductive target with a bulge.

373

points labeled ‘distortion corrected’ are thus produced. To test the success of the final procedure, a model with a flat surface and containing only the crustal conductor was used to calculate the apparent resistivity curve labeled ‘bulge only’ in Fig. 7.

~

1111111

ti~ ~

I

III

IIII~

I

w

I

11111

1~.4,jI IlU

MEXICO SALAD4 BASIN

TE (DC)

LAGUNA

00

I

111115

—

I

i

s..~

>-

It is obvious from the agreement in Fig. 7 between the ‘bulge only’ curve and the distortioncorrected points that tensor stripping works we!! in this case. The overall level of the apparent resistivity values near 50 ~2m is due to the entire 1042-m conductive half-space, from 7.5-km depth to infinity (Fig. 6). The apparent resistivity effect of the bulge itself is a subtle decrease at 10-s period (— 15%). Even so, both the proper and the distortion-corrected results show this apparent resistivity decrease centered at 5-km distance. The two results differ by a maximum of 6% (3 parts in 50). To illustrate the use of the Ray!eigh—FFT in distortion correction of actual MT field data we present a result from northern Baja California, Mexico. The observed TE and TM sounding curves are shown in Fig. 8 from site 10 in the Laguna Salada Basin of the Salton Trough, a region of active crustal spreading. At the longest periods, after 400 s, the TE curve appears to ‘turn over’, which could be indicative of a deep crusta! conductor. However, the conductive basin with resis-

H

10 ~__—

~ ~ ~ ~

SITE 10

I

0

1111111

0.02

0.1

I

-.

_____

.

I

1111111

I

_____

111111

10

I

11111

_____

~III

00

I0002000

< PERIOD (SEC) Fig. 8. Magnetotelluric apparent resistivity period sounding curves at site 10 in the Laguna Salada versus Basin, Mexico and distortion-corrected TE (DC) curve.

—

I 250T

I

I

I

I

I

TM

10 SEC

I x

200

HILL & BASIN

model (2—3000 ~2m) are (self-consistency probably responsible for insufficient convergence checks)

~

0

150

\

~

z

ULGE/ BULGE ONLY~

~

DISTORTION

-~~~~REcTEDX

Ui

a. a.

in the TM results. However, the tensor distortion-

HILL, BASIN ~

—

50

—

4

-

0

I

I

I

1

2

3

I

tivities ranging from 2 to 20 ~2m has significantly affected (or distorted) the deep soundings. The TM results have been distorted by galvanic effects caused by the basin contact with very resistive (— 3000 ~2m) batholithic host rocks and the inductive response of the basin is clear in the TE results. After detailed 2-D gravity modeling along the MT traverse, a basin configuration was derived. The Rayleigh—FFT scheme was then used to MT model the basin and the surrounding mountainous topography (Miele, 1986). This provided the distorted 2-D electric and magnetic fields to implement eqns. (8) and (9). The undistorted model is a 300042-rn flat, homogeneous half-space. The high resistivity contrasts in the geoelectric

I

4 5 6 DISTANCE 1KM)

I

I

I

7

8

9

—

10

Fig. 7. Apparent resistivity versus distance for the TM mode at 10-s period for Fig. 6 models of the bill and basin only; the hill, basin, and bulge; the bulge only; and the distortion-corrected values.

corrected dence thatTE aSalada results deep Basin. conductor in Fig.Removing 8 does provide exist strong beneath evithe Laguna the effect of the basin has produced a very striking descending, long-period branch in the TE sounding curve. This corrected TE curve is nearly parallel to the distorted TM curve (Fig. 8); however, the amount of vertical shift could not have been predicted from this single sounding alone. Three-dimensional effects have not been considered in this field example; they should be evaluated before final interpretations are made.

374

6. Conclusions The Rayleigh— FFT modeling approach can accurately be used to calculate the MT response of many realistic topographic features and subsurface geoelectric structures. The modeling input format is not as tedious as other schemes and by using various self-consistency criteria, valid results can be confirmed and invalid ones rejected. We have used a constant resistivity in this paper to calculate the undistorted electric and magnetic fields. For the correction of some field data we have used a homogeneous half-space resistivity value variable with station location. Magnetotelluric terrain correction using site-dependent homogeneous ground has also been applied by Chouteau and Bouchard (1988). The Rayleigh—FFT algorithm has been successfully used to remove the unwanted effects of topography and near-surface inhomogeneities by distortion tensor stripping of both theoretical and actual field MT results. This paper presents 2D modeling results only. However, an extended version of the computer program (Boersma, 1986) provides a hitherto unavailable code to mode! 3-D Earths with smoothly varying surfaces and interfaces. This is in contrast with the previous block modeling techniques.

Acknowledgments The authors wish to thank John Boersma and Ed Gustafson for their assistance in running the Rayleigh—FFT code on CRAY supercomputers. Personnel at Boeing Computer Services and the San Diego Supercomputer Center were also very helpful. We also thank I.K. Reddy for pointing out a computational error in the first tensor stripping results (Reddig, 1984). The Rayleigh—FFT MT modeling effort has been supported by funds from the National Science Foundation, Grants EAR-8213847 and EAR-8609292. The field result was collected through joint effort with Centro de Investigación Cientifica y de Educación Superior de Ensenada, Mexico. The initial version of this paper was written while one of us (G.R.J.) was a senior Fulbright Scholar, under the auspices of the

Australian—American Educational Foundation, at Macquarie University, Sydney, Australia. We also wish to thank Art Raiche, Ben Sternberg, and Phil Wannamaker whose reviews of the initial manuscript contributed to several improvements in the final version.

References Aki, K. and Lamer, K.L., 1970. Surface motion of a layered medium having an irregular interface due to incident plane SH waves. J. Geophys. Res., 75: 933—954. Andrieux, P. and Wightman, W.E., 1984. The so-called static corrections in magnetotelluric measurements. Presented at 54th Annual Society of Exploration Geophysicists Meeting, Atlanta, Expanded Abstracts, pp. 43—44. Berdichevsky, M.N., Vanyan, L.L., Kuz.netsov, V.A., Levadny, VT., Mandelbaum, M.M., Nechaeva, G.P., Okulessky, B.A., Shilovsky, P.P. and Shpak, I.P., 1980. Geoelectrical model of the J.R., Baikal1986. region. Phys. Earth 22: 1—11. Boersma, Application of Planet. RayleighInter., Scattering Theory to Three-Dimensional Magnetotelluric Modeling. M.S. Thesis, San Diego State University, San Diego, CA. Bostick Jr., F.X., 1986. Electromagnetic array profiling. Presented at the 56th Annual Society of Exploration Geophysicists Meeting, Houston, Expanded Abstracts, pp. 60—61. Chouteau, M. and Bouchard, K., 1988. Two-dimensional terrain correction in magnetotelluric surveys. Geophysics, 6: 854-862. Jiracek, G.R., 1972. Geophysical Studies of Electromagnetic Scattering from Rough Surfaces and from Irregularly Layered Structures. Ph.D. Thesis, University of California, Berkeley, CA. Jiracek, G.R., 1973. Numerical comparisons of a modified Rayleigh approach with other rough surface scattering solutions. IEEE Trans. Antennas Propagat., AP-21: 393—396. Jiracek, G.R., Rodi, W.L. and Vanyan, L.L., 1987. Implications of magnetotelluric modeling for the deep crustal environment in the Rio Grande rift. Phys. Earth Planet. Inter., 45: 179—192. Jones, A.G., 1988. Static shift of magnetotelluric data and its removal in a sedimentary basin environment. Geophysics, 7: 967—978. Kaikkonen, P., Vanyan, L.L., Martanus, ER. and Okulessky, B.A., 1985. Contribution of the surficial effects on the low-frequency magnetotelluric anomaly at the Rheingraben area. Phys. Earth Planet. Inter., 37: 223—227. Kojima, R.K., 1985. Application of the Rayleigh—FFT Technique to Magnetotelluric Modeling. M.S. Thesis, San Diego State University, San Diego, CA. Larsen, J.C., 1977. Removal of local surface conductivity effects from low frequency mantle response curves. Acta Geodaet. Geophys. Mont., 12: 183—186.

375 Lippmann, B.A., 1953. Note on the theory of gratings. J. Opt. Soc. Am., 43: 408. Miele, M.J., 1986. A Magnetotelluric Profiling and Geophysical Investigation of the Laguna Salada Basin, Baja California. M.S. Thesis, San Diego State University, San Diego, CA. Mozley, E.C., 1982. An Investigation of the Conductivity Distribution in the Vicinity of a Cascade Volcano. Ph.D. Thesis, University of California, Berkeley, LBL-15671. Park, S.K., 1985. Distortion of magnetotelluric sounding curves by three-dimensional structures. Geophysics, 50: 785—797. Park, S.K., Orange, A.S. and Madden, T.R., 1983. Effects of three-dimensional structure on magnetotelluric sounding curves. Geophysics, 48: 1402—1405. Rayleigh, Lord (J.W. Strutt), 1896. The Theory of Sound, Vol. 2, 2nd edition, section 272. Macmillan, London. Reddig, R.P., 1984. Application of the Rayleigh—FFT Technique to Topographic Corrections in Magnetotellurics. M.S. Thesis, San Diego State University, San Diego, CA. Reddig, R.P. and Jiracek, G.R., 1984. Topographic modeling and correction in magnetotellurics. Presented at 54th Annual Soc. Expl. Geophys. Meeting, Atlanta, Expanded Abstracts, pp. 44—47. Rokityansky, 1.1., 1982. Geoelectromagnetic Investigation of the Earth’s Crust and Mantle. Springer, Berlin, 381 pp. Shoemaker, C.L., Shoham, Y. and Hockey, R.L., 1986. Interpretation of natural source electromagnetic array data. Presented at the 56th Annual Soc. Expl. Geophys. Meeting, Houston, Expanded Abstracts, pp. 63—65.

Sternberg, B.K., Washbume, J.C. and Anderson, R.G., 1985. Investigation of MT static shift correction methods. Presented at 55th Annual Soc. Expl. Geophys. Meeting, Washington, DC, Expanded Abstracts, pp. 264—267. Stratton, J.A., 1941. Electromagnetic Theory. McGraw-Hill, New York, 615 pp. Vanyan, L.L., 1981. Deep geoelectrical models: geological and electromagnetic principles. Phys. Earth Planet. Inter., 25: 273—279. Wannamaker, P.E., 1983. Resistivity structure of the Northern Basin and Range. In: The role of heat in the development of energy and mineral resources in the Northern Basin and Range Province. Geothermal Resources Council, Special Rep., 13: 343—362. Wannamaker, P.E., Hohmann, G.W. and Ward, S.H., 1984. Magnetotelluric responses of three-dimensional bodies in layered earths. Geophysics, 49: 1517—1533. Wannamaker, P.E., Stodt, J.A. and Rijo, L., 1986. Two-dimensional topographic responses in magnetotellurics modeled using finite elements. Geophysics, 51: 2131—2144. Ward, S.H., 1967. Electromagnetic theory for geophysical application. In: Mining Geophysics, Vol. 2. Society of Exploration Geophysicists, Tulsa, pp. 10—196. Word, D.R., Goss, R. and Chambers, D.M., 1986. An EMAP case study. Presented at the 56th Annual Soc. Expl. Geophys. Meeting, Houston, Expanded Abstracts, pp. 61—63.