RAYINVR v2.0 is an enhanced version of RAYINVR which contains the option to allow forward-modelled rays to be traced through a two-dimensional model with a radial parameterisation. This allows the user to take the sphericity of the Earth into consideration when inverting travel times.
See also: Gorman, A.R. 2002. Ray-theoretical seismic traveltime inversion - Modifications for a two-dimensional radially parametrized Earth. Geophysical Journal International 151: 511-516.
Modifications by: Andrew R. Gorman
Presently at: Department of Geology
University of Otago
Dunedin 9015, New Zealand
e-mail: andrew.gorman@stonebow.otago.ac.nz
Original program written by: Colin A. Zelt
Presently at:
Rice University
Houston, Texas 77005
USA
e-mail: czelt@rice.edu
Controlled-source seismic data are most often modelled with the assumption that we are living on a flat Earth. It is assumed that the "depth" co-ordinate, z, is straight down and that the "horizontal" co-ordinates, x, lie in a straight line along the surface. This approximation is fine for most modelling situations. However, when a seismic survey becomes large enough (e.g., the Russian Peaceful Nuclear Explosions data sets or the Deep Probe data set in western North America), the Earth's sphericity can become quite important. When this is the case, it is necessary to consider that the values of x and z really are referring to radial depth and angular distance along the surface of a sphere.
If the Earth's sphericity is not taken into consideration, then error will be present in the velocity structure of the forward model. This error will increase with depth. In the past, this error has been considered through the use of an Earth flattening approximation (e.g., Müller, 1971) which shifts measured traveltimes through the Earth by an amount equal to the difference between modelled traveltimes though one-dimensional spherical and flat models.
The main concern with this approach is that the approximation does not take lateral variations in velocity into consideration. A better approach is to perform the ray tracing in a radial co-ordinate system.
First, consider how the ray tracing equations are formulated in a `Flat' Earth, i.e., where geometries are expressed using the Cartesian co-ordinate system.
Figure 1: Cartesian GeometryIn this system, the ray tracing equations are a pair of ordinary differential equations (Cervený et al., 1977) (Table 1) cast in terms of `depth', z, and `offset distance', x. For computational efficiency, the orientation of the ray can be taken into consideration to produce two equivalent pairs of ODE's: one pair when the ray is mostly vertical (incidence angle between 45° and 90°), and another when the ray is mostly horizontal (incidence angle between 0° and 45°.)
Travel times are then simply calculated by integrating slowness (reciprocal of velocity) along the raypath, s.
Now consider a more accurate parameterisation of the Earth in cross section: a two-dimensional Earth with radial geometry. In this system, `depth' is converted to radial distance (distance from the centre of the Earth) and `offset distance' between the source and receiver is recast as an angle. Table 2 summarises the equivalents in each parameterisation.
Figure 2 shows the geometry of this new parameterisation.
A new set of ordinary differential equations (Table 3) can now be produced using the equivalents determined in Table 2.
It is not convenient or practical in controlled-source seismology to parameterise an input model in terms of radial distance and angle. However, it is simple to modify the radial equations in order to use the same input parameters (x and z) as for a flat Earth. This is done by including the radius of the Earth as a constant. Table 4 summarises the equivalents of these two parameterisations.
Figure 3 shows the radial Earth parameterisation in terms of x and z.
Finally, two pairs of Radial Earth Ordinary Differential Equations (Table 5) are produced which are parameterised in terms of depth and offset distance.
Note the similarity between these ODE's and the ODE's used for a flat Earth. The only difference is the application of a scaling factor involving the ratio of the Earth's radius to the radial distance.
One new switch and one new parameter have been added to the ray tracing parameters (TRAPAR namelist). These are ispher and rearth.
rearth - radius of the Earth. (default: 6371) This could be modified for equatorial or polar regions where 6371 km is an unsatisfactory approximation for the radius of the Earth. It could also be modified for use on another planet.
In the file rayinvr.com a new commonblock /blksphere/ has been added to include the parameters ispher and rearth.
The RAYINVR subroutine, TRACE, located in file trc.f controls the tracing of a single ray through the velocity model given an initial set of conditions. This subroutine traces the ray in x, z space. In rewriting the code to accommodate radial Earth considerations, the parameterisation of the input model was left the same, i.e., in x and z. However, the radial nature of the Earth means that units of x now get smaller with depth. This must be considered for all RAYINVR routines which use z and x to calculate length or time. The subroutine TRACE has been modified as follows.
TRACE uses one of two Runga-Kutta algorithms to solve the pairs of ordinary differential equations (ODE's) which make up the ray-tracing system. One of these algorithms has error control (subroutine RNGKTA) and the other does not (subroutine RKDUMB). Which one is used is specified in the TRAPAR namelist with the switch ifast. Each Runga-Kutta routine calls an addition subroutine containing the ODE's that it solves numerically. New ODE subroutines have been added (as described above) for radial geometries.
Four new ODE subroutines were written.
RAYINVR contains several supplementary subroutines which assist in ray tracing and in calculating travel times. As mentioned previously, the radial nature of the Earth means that units of x now get smaller with depth. Since the model parameterisation remains in units of x and z then several RAYINVR subroutines which use z and x to calculate length or time must be modified. Most important among these subroutines is TTIME which calculates the travel time from source to receiver for each ray by integrating along the ray path, as described above.
Plotting for the radial Earth case is complicated by the fact that we are dealing with a radial co-ordinate system, but the plotting routines used deal with Cartesian co-ordinates. Several subroutines were developed to transfer values back and forth between the original `x' and `z' parameterisation of radial co-ordinates, and the projection of this co-ordinate system into cartesian co-orinates for plotting to the screen or paper. The orignal parameterisation is generally in units of kilometres, but the transformed values are generally in pixel units.
PLTMOD calls numerous other subroutines to actually perform the plotting functions. Some of these subroutines had to be modified for the radial representation (including SPHRPLT, CURVE, DASHCURVE, and AXISSPHERE.)
The vertical exaggeration of radial co-ordinates takes some careful consideration. The routine is designed to accurately represent the 1:1 curvature of the Earth at any particular radius, using the scale shown for the horizontal axis. The "horizontal" shortening seen with depth is similarly correct.
xin, zin - input co-ordinates in units defined by the plotting routines (pixels)
Cervený, V., I.A. Molotkov and I. Psenvcik, 1977, Ray method in seismology, Univerzita Karlova, Prague.
Müller, G., 1971, Approximate treatment of elastic body waves in media with spherical symmetry, Geophys J R Astr Soc, 23: 435-449.