A FORTRAN 77 program for inverting gravity anomalies of two-dimensional basement structures

Computer* dt Geo*cwmes Vol. 15. No. 7. pp. 1149-1150. 1989 Pnnted in Great B r t t a m All rights reserved

(1098-3004 89 $3.00 -,- 0.00 Copyright ~ 1989 Pergamon Press plc

SHORT NOTE A F O R T R A N 77 P R O G R A M F O R I N V E R T I N G G R A V I T Y A N O M A L I E S OF TWO-DIMENSIONAL BASEMENT STRUCTURES I. V. RADHAKRISHNAMURTHY and S. JAGANNADHARAO Department of Geophysics, Andhra University. Visakhapatnam-530 003. AP., India (Received 12 Juh" 1988; revised I September 1988 ; received lor publicatton 21 March 1989)

INTRODUCTION The usual practice of inverting gravity anomalies of two-dimensional bodies is to replace their cross sections by an n-sided polygon and to determine the locations of the vertices that best explain the observed anomalies (e.g. Johnson. 1969; Enmark. 1981). The initial coordinates of the vertices are assigned and later modified iteratively so as to minimize the differences between the observed and calculated anomalies. The estimation of the initial values is a separate and indeed a critical exercise. This selection determines the convergent solution to the problem. It seems that inversion schemes replacing the twodimensional bodies by a series of juxtaposing prisms, instead of a polygonal cross section, do not require any a priori calculation of the initial values of the parameters that define the outline of the body. This paper presents such an inversion scheme and the related computer program in FORTRAN 77 for determining the density surface such as the basement topography above an assigned depth Z and density contrast a. The method does not require input of initial values of any other parameters. It also is applicable for determining structures with a flat top or a flat bottom.

ly along a profile. It may not be possible to get the anomalies free from a regional component. The program. therefore, determines the regional anomaly as well. The input to the computer consists of gravity anomalies Ag. sampled at equal intervals along a profile, the depth Z, over which the structure or the density surface is to be worked out and the density contrast a. The density surface or structure is equated to a series of juxtaposed prisms (Fig. !) one below each anomaly point, and the depths Z T to their top are determined. The basic approximation used is that the vertical gradient of gravity anomaly of any prism is constant over the depth extent of the prism or the incremental part of it. and that its contribution to the anomaly at any point is a simple multiple of the gradient and thickness or the increment to it. The anomaly profile is assumed to have covered the structure completely, so that the rel'ief below the first and last points becomes zero. The profile is assumed to have a linear regional trend defined by Ax + B. The gravity anomaly Ag(k) = Ag(x,) at any point P(x~) which is the result of the surface or structure is replaced by prisms (Fig. i) and can be written from the gravity anomaly equation of a dike of Rao and Murthy (1978): N- I

Ag(k) =

THE METHOD The program determines depths to the top of the basement surface below each point of gravity anoma-

0

0

•.m--"ClI --e.. 0 0

j~

0

X(i) 0

Z!(I} I

ToDography

~ [F,(z)l~,n,, + Ax, + B 1=2

with

0

Fk(:) = 2Ga {z [arctan((xk + dx/2)/z)

0

PCXk ) 0

O

O

]

ReoLacing the ~opogrophy

tteration

Figure I. Gravity topography and its replacement by prisms. 1149

(1)

1150

Short Note arctan((.,q - dx 2):)]

-

[F

- (.q - dx2)ln(ix~ - dx 2): + ..z)]}

~-~1z z ~ IIZT IH

A~ -

+ 0.5 [(.q + dx 2)In ((x~ 4- dx :2): + .": )

= [ SFA(:---~)] dZT(i) L 8 : J z r ....

(8)

and hence

(2) where dx is the station spacing, ZT(i) is the depth to top of the surface below the ith station and G is the gravitational constant. Notice that uzr., . - r, F ^ ( .-,Jz

[&(z)]zzr,,, =

(3)

A prism lying below the mean depth (ZT(i) > Z ) therefore produces a gravity effect which has an opposite sign to the one produced by a prism lying above the mean depth (ZT(i) < Z). Thus the difference in sign in gravity contributions from the masses lying above and below the mean depth is incorporated automatically. If ZT,(i). i = 2 to N - 1. are the depths to the surface at the end of thejth iteration and A, and B, are the coefficients of the regional anomaly, then the anomaly of such a model is given by

dg(k) =

'7--":

L

c:

dZT(i) + d.4 x~ + dB. Jzr.~ (9)

In Equation (9). xA is a known parameter and I?Fa(:}?Zlz~,~ is calculated by Equation (5). dg(k) is the difference between the observed and calculated anomalies. Equations of this type are framed, one for each anomaly point, and solved for values ofdZT(i). dA and dB through minimizing the objective function 'r~= t dg(k)-" by employing Marquardt's (1963) optimization method. The normal equations are written as;

- -

( I + (2)a,

dg(k). ~ .

,,.~ ..... (10)

From Equation (2). H"A(:) ~Tz

=

where the constant 6 = I for i = j, and 0 otherwise. :. is the Marquardt's damping factor.

-- 2G~ [arctan[(.xA + dx/2)/:]

~t

- a r c t a n [ ( x x - dx/2)/:]].

tl,,,

This equation resembles the gravity anomaly of a thin plate model, the only difference being the negative sign. The present scheme implicitly consists of two phases, one for initialization and the other for refinement. For initialization, it is assumed that each prism will contribute to the anomaly a component equal to the product of its thickness (Z - ZTt (i)) and the rate of change of anomaly with Z T = Z. Thus Ag(k) =

'Z' [dFx(z) 3 ~ (7. - ZT,(i)) ,-: L

c:

1

=

I to,V - 2

A

and a~ = B so that kdg(k) tTdg(k)

for/' =

F?Fx (:)]

-- x(k)

ItoN-2

and

~gl ,¢ _ I

?dg(k) (6)

=

1.0.

~(-/3,

Jz

+ A,xx + B,. from which the initial values, ZT,(i), A,, and B,, of the depths to the surface and the coefficients of regional anomaly are calculated. The difference dg(k) between the observed anomaly and the one calculated at the end ofjth iteration is written from Equations (I) and (4) as dg(k) =

= ZT(i + I) f o r i =

(5)

E [Fx(:)]zr:,~ zr~,~ + dA xx + dB,

(7)

1-2

where dA and dB are the increments to be given to A, and B~. Ifthe increments dZT(i) = ZT(O - ZTj(i) to be given to the calculated depths are small, we can ~rite that

The iteration is terminated when the Marquardt's damping factor 2 becomes large and exceeds a defined value say 15.0, the objective function is reduced to an allowable error, or the specified number of iterations are completed. It may be observed that even the exercise of initialization actually becomes a built-in part of the iterative scheme, by initially setting ZT(i) = Z. i = I to N, and A = 0 and B --- 0. The anomaly calculated by such a model is zero. Therefore the errors initially are set to the observed anomalies resulting in the equality of Equations (9) and (6). It may be construed that either the first iteration in the suggested scheme effectively calculates the initial values of ZT. A, and B, or that the method starts by equating the initial relief to zero,

Short Note

1[ 51

Observed

"~ "~ IC

'~

Cokculo~ed ~ n o m a ty t 1Okra

~

~

"t~,,,,.

Basement structure Q. O Prlsms representing basement structure

Figure 2. Interpretation of gravity profile over Raguba Field. Sirte Basin. Libya. Gravity profile is after EI-Batroukh and Zentani (1980). In essence, the iteration scheme consists of the following steps: (I) To equate ZT(i) to Z for all values of i = I to N, dg(i) t o A g ( i j f o r i = I t o N , a n d A a n d B t o z e r o , (2) To calcuhttc [,'FA(:)/(":lzn,) for i --- 2 to N - I and t\)r each value o f k = I to N with the help of Equation (5) and s a k e for the increments dZT(i), i = 2 to N - I and d,.I and dB by Equation (9). (3) To add the increments to the respective paralllelers,

(4) To calculate Ag<,i (k) front Equation (4) and then dg(k)= A g ( k ) - Ag<~l(k) for all values of k = I to N a n d , (5) To repeat steps 2, 3, and 4. COMPUTER PROGRAM AND AN EXAMPLE OF INTERPRETATION The computer program coded in F O R T R A N 77 for tracing the two-dimensional topography is given in the Appendix. The input consists of depth above which the density surface is to be traced, number of observations, station interval, observed anomaly and the density contrast. Limiting depths of the surface also may be specified. These limiting depths help convergence to a proper solution and can be selected arbitrarily. When no information is available on the structures, ZTLT. the minimum depth to surface can be equated to zero. whereas ZBLT, the maximum depth can be 3Z or a few kilometers. This program makes use of two function subprograms GPLATE and GPRISM and one subroutine SIMEQ. Subprograms GPRISM and GPLATE calculate gravity anomaly and its vertical derivative of a prism. Subroutine SIMEQ is used to solve a set of simultaneous linear equations. The printout consists of the gravity anomalies and the depths to the interface at the end of each iteration, the printout at the last iteration being the best fitting interpretation. Figure 2 illustrates the published gravity data over

Raguba Field, Sirte Basin, Libya (Et-Batroukh and Zentani, 1980) and its interpretation, with the algorithm presented here, for a two-dimensional basement structure. The depth of the basement and its density contrast are assumed to be 5.7 km and 0.2 gm/ cm ~ respectively. The limiting depths Z T L T and ZBLTare assigned arbitrarily 0 and 10 km respectively. Obviously, the basement depth cannot be < 0. The most appropriate value for ZBLT in the present example is 5.7 kin, but a value of 10kin is used only to demonstrate the flexibility on its selection. The interpretation predicts an unaccounted linear regional anomaly defined by - 0 . 2 x + 1.25mgal. Figure 2 also shows the positions of the interpreted prisms and the corresponding basement structure. The fit between observed and calculated anomalies is excellent. The major basement features as obtained in this interpretation broadly agree with those reported by EIBatroukh and Zentani.

A¢'knnwled¢ment--Th¢ subroutine SIMEQ is the modified version of the subroutine SIMQ of the scientific subroutine package of the IBM available with the computer centre. Andhra University for solution of simultaneous equations. REFERENCES

Enmark, Th., 1981. A versatile interactive computer program for computation and automatic optimization of gravity models: Geoexploration, v. [9, no. I. p. 47-66. Johnson, W. W., 1969. A least squares method of interpreting magnetic anomalies caused by two-dimensional structures: Geophysics, v. 34. no. I. p. 65-74, Marquardt, D. W., 1963. An algorithm for least-squaresestimation of non-linear parameters: Jour. Sac. Indust. Appl. Math., v. II no. 2, p. 431-441. Rao, B. S. R., and Murthy. I. V. R.. 1978, Gravity and magnetic methods of prospecting: Arnold-Heinemann (India) Pvt. Ltd., AB,9 Safdar jang Enclave. New Delhi, 390 p. Salah. L., EI-Batroukh. and Zentani, A. S., 1980. Gravity interpretation of Raguba Field. Sirte Basin, Libya: Geophysics v. 45, no. 7, p. 1153-1163.

1152

Short Note

APPENDIX * PROGRAm'; G R 2 D I N DETERMINES TWO-DIMENSIONAL DENSITY * ABOVE A ~NOh%" DEPTH FROM GK~VITY ~N'OMALIES.

SL~FkCE

* *

INPUT * * * * * * * *

Z

:

DEPTH ABOVE ~ I C H THE STRUCTURE SL~FACE IS CALCUL%TED, k-M

NOBS DX G DC ZTLT

: : : : : :

~MBER OF SAMPLED.~NOMALIES STATION I.';TERVAL, ~

: :

ANOMALYOF STRUCTURE CALCULATED, IN .~JALS * DEPTH TO INTERFACE, K~

ZBLT

ML%SL~ED DENSITY MINIMUM ~Xl~M

ANOMALY, IN MGALS CONYRAST,GM/CC DEPTH TO INTERFACE, DEPTH TO INTERFACE,

KM K~I

OR DENSITY* * * * * * * *

OUTPUT GCAL

•

ZT

IMPLICIT REAL*8 (A-H,0-Z) DIMENSION X(40),G(40),ZT(40),ZTD(40),ERROR(40),DER(40), *p(40,41),PD(40o41),DZ(40),GCAL(40) REAL LAMDA,LAMDAI READ(*,B00)Z,DC,DX,NOBS WRITE(*,[OOO)Z,DC READ(*,BOI)ZTLT,ZBLT

2

READ(*,OO2)(G(I),I-I,NOBS) IF(ZBLT.EQ.O.O)ZBLT-2*Z ABSDC-DABS(DC) DO 2 I-I,NOBS X(t)-(~-l)*DX N~MMI)A-[ CON-I.0 DXBY2"DX/2 NOBSI-NOBS-I NOBS2-NOBS-2 FUNCTIoO,O DO 5 K-I,NOBS ZT(K) " Z

ZTD(K)- Z ERROR(K)-G(K) 5

FUNCTI'FUNCTI+ERROR(K) ~ 2

AER-NOBS*0.000I A -0.0 B

"0.0

WRITE(*,910) FUNCT! DO I00 ITER-I,20 DO 15 M-I,NOBS+I DO 15 I-I,NOBS 15

P(I,M)'O.O DO 40 K'I,NOBS DO I0 I'2,NOBSI

XI-X(K)-X(I) 2TI'ZT(1) I0

DER(I-I)--GPLATE(XI,DXBY2,ZTI,ABSDC)

DER(NOBSI)-X(K) DER(NOBS)-I.O DO 20 I-I,NOBS

20 25 40 42

DO 20 M-I,NOBS p(I,M)-P(I,M)+DER(1)*DER(M) DO 25 I-I,NOBS P(I,NOBS+I)'P(I,NOB S+I)+ERROR(K)*DER(1) CONTINUE DO 45 I'I,NOBS DO 45 K'I.NOBS+I IF(I.EQ.K)THEN

PD(I,K)-p(I,K)*CON ELSE pD(I,K)'P(I,K) 45

ENDIF CONTINUE CALL SIMEQ(PD,DZ,NOBS ,KS) IF(KS.EQ.I)THEN

WRITE(*,IO04) STOP

Short Note

50

60 80

90

95

98

100 800 901 801 802 905 906 907

908 910 911 912

ENDIF DO 50 I'2,NOBSI DDZ'DZ(I-I) ZTD(1)-ZT(1)+DDZ IF(ZTD(I).LT.ZTLT)ZTD(I)'ZTLT IF(ZTD(1).GT.ZBLT)ZTDfl)'ZBLT AD - A+DZ(NOBSI) BD - B+DZ(NOBS) DO 80 K=I,NOBS GCAL(K) - 0.0 DO 60 I-2,NOBSI XI = X(K)-X(I) ZTI -ZTD(1) DT-GPRISM(XI,DXBY2,ZTI,Z,ABSDC) GCAL(K) = GCAL(K)+DT GCAL(K) - GCAL(K) +AD~X(K)+BD ERROR(K)- G(K)-GCAL(K) FUNCT-O.0 DO 90 I=I,NOBS FUNCT=FUNCT+ERROR(1)**2 IF(FUNCT.LE.FUNCTI)THEN DO 95 K-2,NOBS[ ZT(K)-ZTD(K) A-AD B=BD FUNCTI=FUNCT LAMDAI=LAMDA NLAMDA-NLAMDA-I IF(NLAMDA.LE.O)NLAMDA'I LAMDA ° 0.5~(2**(NLAMDA-I)-I) CON " I + LAMDA ELSE IF( LAMDA .GT. 15. O)TIIEN WRITE( ~, lO0[) STOP ENDIF LAMDA l° LAMDA NLAMDA'NLAMDA+l LAMDA " 0.5*(2**(NLAMDA-[)-|) CON " | + LAMDA WRITE(*, 906) WRITE(*,914) ITF~ WRITE(~,912)FUNCT,LAMDAI,LAMDA GO TO 42 ENDIF WRITE(*,906) WRITE(~,908) ITER WRITE(*,907) DO 98 K'[,NOBS WRITE(*,90I) X(K),ZTD(K),G(K),GCAL(K),ERROR(K) WRITE(R,905) WRITE(*,911)AD,BD WRITE(*,913) FUNCT,LAMDAI IF(FUNCT.LE.AER)THEN WRITE(~,I002) STOP ENDIF CONTINUE WRITE(*,[003) FORMAT(3FI0.2,I5) FORMAT(SX,5(FI0.Z,5X)) FORMAT(2FI0.2) FORMAT(8F10.2) FORMAT(5x,;5(IH-))

FORMAT(l~1~) FORMAT(SX,TS('-')/IOX,'DIS-',IOX,'DEPTH TO',6X,'OBSERVED',6X, *'CALCULATED',gX,'ERROR',IIOX,'TANCE',8X,'INTERFACE',6X, *'ANOMALY',9X,'ANOMALY',/5X,75('-')) FORMAT(5X'AT THE END OF ITERATION NO.',I3) FORMAT(SX,'OBJECTIVE FUNCTION'',FI0.2) FORMAT(5X,'A-',3X,FT.2/ 5X,'B-',3X,FT.2) FORMAT(SX,'THE OBJECTIVE FUNCTION I S ' , F I O . 0 , 1 * 5X,'THE CURRENT VALUE OF IAMDA IS ',FI0.31 5X,'AS THERE IS AN INCREASE IN OBJECTIVE'/ * 5X,'FUNCTION LAMDA IS INCREASED TO',FI0.3,1 * 5X,'AND CALCULATION IS REPEATED')

1153

1154

$ h o n Note 913

FORMAT(SX, 'OBJECTIVE FUNCTION=' ,FI0.2/ 5X, ' LAMDA-', 13X,F10.2) FORMAT(SX,' ITERATION N O . ' ' , I 3 / 5 X , 1 3 ( ' - ' ) ) FORMAT(SX,'I.~VERSION OF GRAVITY ANOMALIES OF A TWO-dIMENSIONAL',/, *SX, ' STRUCTURE AT A KNOWN DEPTH' , F 5 . 1 , ' KH. ' , / , 5X, ' THE ASSUMED DENSI tTY CONTRAST I S ' , F 5 . 2 , 'GM/CC.' , / / / ) FORMAT(///,5X,'LAMDA IS INCREASED TO AN UNUSUALLY LARGE VALUE. BET ~*TER'/,5X, 'SOLUTION THAN THAT IN THE LAST PRINTED TABLE MAY NOT'/SX t, 'BE FOUND. ALSO CHECK FOR POSSIBLE ERRORS IN INPUT DATA. ') FORMAT(///5X, 'THE OBJECTIVE FUNCTION IS LESS THAN THE ALLOWABLE ER tROR. '/5X, 'THE LAST PRINTED TABLE PROVIDES DEPTHS TO THE DENSITY S tURFACE OR STRUCTURE. ') FORMAT(///SX, 'THE INTERPRETATION MAY NOT IMPROVE BY PERFORMING FUR *THER ITERATIONS. '/SX,'THE LAST PRINTED TABLE PROVIDES THE DEPTHS T tO THE DENSITY SURFACE OR STRUCTURE. ' ) FORMAT(///5X,'ILL-CONDIT10NED MATRIX. CHECK DATA. ' / ' THE PROGRAM tIS TERMINATED. ' ) STOP END t

914 1000

1001

1002

1003

1004

It e t t t t e t

***eteett*eetttttt

t e fttemetttmet

tettttttttm~t

* FUNCTION SUBPROGRAM GPLATE CALCULATES GRAVITY * ANOMALY OF A THIN PLATE MODEL.

*

INPUT

t * t

X T Z DC

: : : :

DISTANCE TO POINT OF CALCULATION * HALF WIDTH t DEPTH DENSITY CONTRAST,GM/CC *

~ t * * ~ t e t ~@ectet e * e t e t e e Jt* ~ e c ~ e t t t t Jtt f e t e t ~ t e t t

2 3

4 5

* *

e~**~,~

FUNCTION GPLATE(X, T, Z, DC) IMPLICIT REALt8(A-H ,O-Z) TWOGS'I3. 333333"DC XPT-X+T XMT-X-T IF(Z,NE.0,0)THEN G PLATE -DATAN (X PT/Z )- DATAN (XMT/Z ) ELSE IP(DABS(X)-T)2,3,4 GPLATE-3. 14159265 GO TO 5 GPLATE'I .570796325 GO TO 5 G PLATE-O. 0 ENDIF GPLATE'GPLATEtTWOGS RETURN END

tiff ~ t t W i t f f t @ # i t t t e ~ # # W t ~ t t e ~ t t ~ t t # ~ * t t t ~ t t t * e ~ t t

~tttt~tet~t

~

, FUNCTION SUBPROGRAM GPRISM CALCULATES GRAVITY ANOMALY ANOMALY OF A VERTICAL PRISM/DYKE.

*

*

*

*

INPUT

t

.....

t t * t

X T ZT ZB DC

*

" DISTANCE TO POINT OF CALCULATION FROM THE * ORIGIN LYING ABOVE THE CENTRE OF THE PRISM.* t HALF WIDTH OF THE PRISM. * t DEPTH TO TOP OF THE PRISM. • DEPTH TO BOTTOM OF THE PRISM. : DENSITY C0N~T,GM/CC. t

ff~tttCrt~t~att~tl~tCtt#~Att~t~tt~t~ett

FUNCTION GPRISM(X,T, ZTOP, ZBOT,DC) I M P L I C I T REAL*8 (A-H, O-Z) TWOCS-13.333333"DC IF(ZTOP. EQ. ZBOT)THEN GPRISM-O. 0 ELSE IF(ZTOP. LT. ZBOT) THEN ZT'ZTOP ZB'ZBOT C'TNOGS ELSE ZT=ZBOT ZB-ZTOP C--TWOGS ENDIF XPT-X+T XHT-X-T

t~t t~tffttt~@~tt*~te@

Short Note

2 3 4 5

I155

IF (ZT. NE. 0.0)TH~N AI =DATAN (XPT/ZT)-DAT~¢(XMTI ZT) ELSE IF(DABS(X)-T)2,3,4 AI=3.14159265 GO TO 5 AI=I.579796325 GO TO 5 AI'0.O ENDIF A2=DATAN(XPTIZB)-DATAN(~MTIZB) TEP~HI-ZB*A2-ZTtAI RI=XPT**Z+ZT**2 IF(RI.EQ.O.0)RI=0.1E-05

R2"XMT**2+ZT**2 IF(R2.EQ.0.O)R2"0.1E-05 R3"XPT**2+ZB**2 R4-XMT*~2+ZB**2 TERM2=XPT~DLOG(R3/RI)-XMT*DLOG(R41R2) GPRISM= (TERM I+0.5tTERM2) ~C ENDIF RETURN END * SIMEQ IS MODIFIED VERSION OF THE IBM SUBROUTINE * SIMQ FOR SOLUTION OF N SIMILTANEOU5 EQUATIONS. P IS THE (N+I)*N MATRIX CONTAIING COEFFICIENTS OF EQUATIONS. THE VALUES OF PARAMETERS ARE STORED IN * B AND TRANSFERRED TO THE MAIN PROGRAM.

3 4

28 31 35 36

56

55

61

SUBROUTINE SIMEQ(P,B,N,KS) IMPLICIT REAL~f(A-H,O-Z) DIMENSION P(40,41),B(40),A(1600) I'N+I DO 3 II'I,N DO 3 12=I,N 13-(II-I)*N+I2 A(I3)'P(12,11) DO 4 14=I,N B(14)-P(14,1) TOL=O.0 KS-0 JJ--N DO 68 J=l,N Jy-j+! JJ=JJ+N+l BIGA-O IT-JJ-J DO 311=J,N IJ=IT+I IF(ABS(BIGA)-ABH(A(IJ))) 28,31,31 BIGA-A(IJ% IMAX=I CONTINUE IF(ABS(BIGA)-TOL) 35,35,36 KS=I RETURN II'J÷N*(J-2) IT=IMAX-J DO 56 K=J,N II'II+N 12"If+IT SAVEoA(II) A(II)'A(I2) A(I2)-SAVE A(II)-A(II)IBIGA SAVE-B(IMAX) S(IPt~J()'S(J) B(J)-SAVE/BIGA IF(J-N)55,71,55 IQS=N*(J-I) DO 68 IX-JY,N IXJ=IQS+IX IT=J-IX DO 61JX'JY,N IXJX'N*(JX-I)+IX JJX-IXJX+IT A(IXJX)=A(IXJX)-(A(IXJ)*A(JJX))

* * t *

Short Note

1156 68 71

B(IX)-B(IX)-(B(J)*A(IXJ)) N~t=N-I IT'N*N DO 81 J=I,NY IA=IT-J IB=N-J IC=N DO 81 K'i,J

B(IB)=B(IB)-A(IA)*B(IC) 81

IA=.~A-N IC=IC-I RETURN END