aniprop.f -- a standalone program to generate surface waves in layered anisotropic structures

click title to see Google Drive directory JParkCodes

Sportsfans,

The erstwhile version of aniprop.f that was (is?) available on my website has been superceded by an all-in-one version that computes dispersion curves and synthetic seismograms all in one program. Adapted and debugged on an Intel Mac in gnu fortran.

 aniprop_072307.f  in Google Drive directory JParkCodes

(after a 24 March 2008 bugfix is now aniprop_032408.f)

The code computes surface wave modes at evenly-spaced frequency points between 0 and 0.5 Hz. You can change this by tinkering with the code. It computes synthetics for a source at a specified downrange distance and depth, and plots them using PLOTIT. The dispersion values are output to files out_cvel and out_gvel (phase and group velocities, natch), and the synthetics are written out to SAC-format files.

c aniprop - program to calculate propagating modes of anisotropic layer
c writes a file of dispersion curves at evenly-space freq points
c
c reads a layered model like the following file (ignore the "c "s)
c
c K&H SoCal model, deep crust horizontal anisotropy TITLE
c 3 # OF LAYERS OVER HSPACE
c 45 45 THETA,PHI ORIENTATION ANGLES
c 4000 5500 0.06 0.00 3175 0.03 2600 FOR SYMMETRY AXIS
c 0 0
c 27400 6300 0.00 0.00 3637 0.00 2800 DEPTH (M), VP (M/S), "B", "C",
c 90 45 VS (M/S), "E"
c 32400 6800 0.04 0.00 3925 0.02 2900 B,C,E ARE ANISOTROPIC PARAMETERS
c 0 0
c 60000 7800 0.00 0.00 4500 0.00 3200 NOTE: HSPACE MUST BE ISOTROPIC
c

One quirk of this code is misbehavior for isotropic models. The code solves for Rayleigh and Love dispersion simultaneously, because the two polarizations are hybridized in anisotropic models. Paradox: for a perfectly isotropic model, the rootfinder sometimes misses roots, typically at occasions where the Rayleigh and Love dispersion curves cross. The bookkeeping for solutions is not very rigorous, so loss of a mode can happen. With anisotropy, the hybridization of Rayleigh/Love motion seems to repel the roots and keep some minimal separation. As far as I can tell, the total effect on actual synthetic seismograms is negligible.

Wave propagation in a truly 1-D lossless model leads to some interesting reverbarations that one rarely sees in nature. With all the overtones included, you can see a succession of bouncing S waves that precede the main surface wave pulse. In Earth's crust, such phases would scatter away their high frequency energy from surface topography and sedimentary layers. There is a loop for mode-summing that can be altered to eliminate the overtones, leaving the fundamental Love and Rayleigh waves: (lines 681-3)

c loop over overtones at a certain frequency
do iov=1,maxbr
c do iov=1,2


Notes on the angle conventions for w-hat, the axis of symmetry:

In the anisotropic reflectivity code, subroutine zzget *assumes* a coordinate system in which z is down (anti-vertical), x is the radial direction, and y is anti-transverse. Therefore, the position angles theta,phi for w-hat are tilt relative to down, and azimuth defined as a rotation from x towards y. This rotation is CCW if viewed from below, and CW if viewed from above. Since w-hat and -(w-hat) define the same axis of symmetry, the position angles *also* can be defined as theta=(tilt from vertical) and phi=(rotation from anti-x (anti-radial) towards anti-y (transverse)). Viewed from above, this phi rotation is CW, and defines the strike of w-hat relative to the arrival azimuth of the wave.

In order to compute seismograms for a variety of back-azimuths, the synthetic code genrecd_az.f accepts a layered model with w-hat position angles defined as theta=(tilt from vertical) phi=(strike CW from N). For an event at back-azimuth psi (CW from N), the code rotates w-hat from geographic coordinates to ray-based coordinates before passing it to subroutine zzget. If a wave arrives at back-azimuth psi, the strike of the axis of symmetry w-hat relative to its arrival azimuth is phi'=phi-psi. The code performs this rotation with this code in subroutine matget, for w-hat azimuth "az":

caz=dcosd(az)

saz=-dsind(az) ! sin(-az)

do n=1,nlp

ww(3)=w(3,n)

ww(1)=w(1,n)*caz-w(2,n)*saz

ww(2)=w(1,n)*saz+w(2,n)*caz

...

In this manner, the axes of symmetry of the model, saved in array w(.,.), are never modified.

Journal Reference:

Park, J., Surface waves in layered anisotropic structures, Geophys. J. Int. v126, 173-184, 1996.