Sven Ivansson Feb. 23, 1998
National Defence Research Establishment (FOA)
S-17290 Stockholm, Sweden
email: sveni@sto.foa.se
General introduction to the program rpressfw
--------------------------------------------
The program rpressfw belongs to the rpress package. It can be used for cw
("continuous wave" with a monofrequency time dependence) computations as
well as computation of transients. The underlying theory is described in
the references. A range-independent (laterally homogeneous) fluid-solid
medium is assumed, with a fluid region interspaced between two (possibly
vanishing) solid regions. The lower and upper boundaries may independently
be either free, rigid, or of "homogeneous half-space" type (either fluid or
solid). In addition, an entirely solid medium (without fluid) can be
handled. The wave-field is created by a vertical array of symmetric point
sources. These point sources can be located in the fluid as well as solid
regions, as can the receivers. An isolated mode could be excited by
implementing the theoretical result of Ref. 5. The technique of Appendix B
in Ref. 6 to improve computation speed in the presence of a source array by
using local particular solutions has not been implemented in rpressfw. (It
has been implemented in a program called rpressmulti, however.)
Units are s for time, km for length and "10**12 kg" for mass. However, dB
is "re 1 m" in the way explained below. The output concerns the particular
field component as specified by the parameter KFWTYP, see below. According
to the chosen units, displacement components will be given in "km" and
stress components will be given in "(10**15 N)/(km**2)".
For the cw computations, an underlying harmonic time dependence s(t)
according to "s(t)=exp(-i*omeg*t)" is understood, omeg>0 being the angular
frequency and t being time. For computation of transients, the source
wavelet s(t) can be arbitrarily specified; the same wavelet is assumed for
all point sources in the source array. In any case, a point source of
strength unity (AMSOURC() below equal to one) has a diagonal moment tensor
with all three diagonal elements equal to " s(t) " "(10**15 N)*(km)". The
output for the cw computations concerns the complex field, and phases are
given in addition to the amplitudes.
For typical applications, point sources of strength unity would be
unrealistically strong. Hence, weaker sources should be specified or one
would have to scale down the results. When dB output is requested, it might
be useful to adapt the source strengths to the particular convention for dB
output described below.
As specified in the information file read.me in this directory rpress, there
are other programs in the rpress package. Some important general remarks on
the implementation are also given there.
Input to the program rpressfw
-----------------------------
The only formal difference from rpress is that KFWTYP and IYSTAB are added
at the beginning.
KFWTYP
IYSTAB
input starting with LX LI(TOP) ICWORKDIM, see the information file rpress.txt
Explanation of the parameters
-----------------------------
KFWTYP specifies the field component that is to be computed for a
depth-range (z-r) coordinate system with depth z increasing downwards.
There are four options: 1 for sign-changed horizontal displacement, 2 for
sign-changed vertical displacement, 3 for sign-changed tangential stress in
horizontal cuts, 4 for sign-changed normal stress in horizontal cuts. Please
note that the positive z direction is downwards and that the field components
will be sign-changed! When KFWTYP=4, the output in the fluid for a point
source of strength unity at the matching depth will be consistent with the
pressure output from rpress in this way. The horizontal displacement is left
uncomputed within the fluid and set to zero there.
IYSTAB concerns stabilization during the backward propagation. Stabilization
will only be performed if IYSTAB is nonzero. A negative value will imply
ordinary stabilization (section 5 of Ref. 6), whereas a positive value will
imply modified stabilization (Appendix A of Ref. 6).
For an explanation of the remaining parameters, see the information file
rpress.txt in this directory rpress. We only need to mention a few things:
* In the rpressfw case, the reference pressure "pm0" is as if rpress were
run with a point source (of strength unity) at our matching depth. For
IDB>=0, all cw field results will be given relative "pm01m=1000*pm0",
displacement as well as stress components. This convention should be
carefully observed.
* The option with a nonzero INTREST has not been implemented.
* Concerning the Bessel/Hankel functions specified by IBESSARR() and
IBESSTAIL(), rpressfw may involve J1 instead of J0 and H1(1) instead of
H0(1). This happens precisely when KFWTYP equals 1 or 3. (J0 and H0(1) are
still used when KFWTYP equals 2 or 4.)
Format for the medium infile FILM
---------------------------------
ZSOURC0 ZSOURC1 NS NSHC
NSOURCFW
do i = 1, nsourcfw
NSOFW0(i) AMSOURCFW(i)
end do
AA BB RRHO DBLA DBLB
a number of data blocks of the following type
ZZ IIFSMTYP IIMHQANI
AA BB RRHO DBLA DBLB NSH NH
Explanation of the format for the medium infile FILM
----------------------------------------------------
The vertical source array consists of NSOURCFW individual symmetric point
sources. Point source number n is located at the bottom of layer NSOFW0(n),
referring to the layer numbering prior to possible splitting as governed by
the parameters MHQ0,FREQMHQ0 and the outer extrapolation. The NSOFW0() must
be positive integers and the medium parameters must be continuous at the
corresponding depths. The relative strength (and phase) of point source
number n is given by the complex number AMSOURCFW(n). Nonreal values of
AMSOURCFW() should not be used for calculation of transients. In the
transient case, it seems appropriate to introduce time shifts (rather than
phase shifts) among the point sources, but this has not been implemented.
(There is an option of this kind in a program called rpressmulti, however.)
It follows that the outgoing symmetric displacement field in the vicinity
of point source number n appears as
AMSOURCFW(n) * exp(i*omeg*R/alpha) *
( R**(-2) - i*omeg/(alpha*R) ) / (4*pi*rho*(alpha**2))
where omeg is the angular frequency and R is the spherical distance to the
source in question, while alpha and rho are the medium P-velocity and
density parameters at the source. As mentioned, the positive z direction is
downwards. Please note, however, that the field components will be
sign-changed in the output! (See the previous explanation of the parameter
KFWTYP.)
In other respects, things are essentially as for the medium infile FILM to
rpress, see the information file rpress.txt in this directory rpress. We
only need to make a few remarks:
* Unless specified among the NSOURCFW point sources, there will not be a
source at the matching depth.
* The receiver depths need not be restricted to the (interior) fluid
region, they may appear in the solid as well.
* There is an additional option to handle an entirely solid medium (without
interior fluid). This option is specified by NS=NSHC=0. In this case, one
must take ZSOURC0=ZSOURC1="the initial ZZ", and the matching depth will
be at this depth.
Format for the outfiles
-----------------------
The format is the same as for the corresponding outfiles from rpress, see
the information file rpress.txt in this directory rpress.
Subroutine directories
----------------------
The source code for the program and its subroutines can be found in the
following directories:
rpress/cpm (program and main routines)
rpress/cpm/single (routines that have different "multi" versions)
rpress/cpm/fw ("full-wave" routines)
rpress/int (general integration and extrapolation routines)
rpress/rut (general routines)
rpress/bess (Bessel- and Hankel-function routines)
rpress/slatec (slatec routines for complex Airy functions)
rpress/ik (routines to allocate space for arrays)
Example
-------
We give a simple example of a VSP (vertical seismic profiling) computation
for a fluid-solid medium with three regions, solid-fluid-solid. This example
does not appear in any of the references, and we discuss it in general terms
before providing the corresponding input to rmodpfw.
At the surface (depth z = 0.000 km) there is a 0.001 km thick ice coverage,
and below the ice there is 0.024 km of water. Below the water there is a
layered solid half-space. The medium parameters are specified in a table:
depth z P-veloc S-veloc rho absorption_P absorption_S
ICE 0.000-0.001 3.500 1.600 0.98 0.1 0.10
WATER 0.001-0.025 1.500 - 1.00 0.0 -
BOTTOM 0.025-0.049 2.000 1.200 2.20 0.075 0.15
0.049-0.061 4.000 2.300 2.65 0.60 0.12
0.061-0.075 3.000 1.600 2.40 0.67 0.13
0.075- 6.500 4.000 3.20 0.0505 0.10
Table: Medium parameters for the solid-fluid-solid medium of the VSP example.
Depths are given in km, velocities in km/s, density in (10**12 kg / km**3) =
(kg/dm**3), and absorption in dB/wavelength. The layers are homogeneous. At
the surface there is a free boundary, whereas the bottom below the depth z =
0.075 km is assumed to be a homogeneous half-space.
In a vertical borehole at the bottom, there are 30 vertical geophones with a
spacing of 0.002 km. The uppermost and lowermost geophones are at the depths
z = 0.027 km (i.e. 0.002 km below the water) and 0.085 km, respectively. A
symmetric point source located 0.001 km below the water (z = 0.026 km) and
at a horizontal distance of 0.010 km from the borehole emits a transient
pulse with dominant frequency of the order of 500 Hz.
The resulting offset VSP section was computed numerically with the compound
matrix method as developed in Ref. 6. For this application, it is of course
essential that the wave-field can be computed at arbitrary depths within a
solid region and not only at its surface. Frequency-wavenumber synthesis
was used with 69 different frequencies covering the frequency interval
0-1000 Hz. For each frequency, adaptive integration over the slowness axis
was used (Ref. 1). In total, a large number of boundary-value problems had
to be solved.
The 30 computed traces are shown in a figure that can be printed from the
postscript file vspfig.ps in this directory rpress. The horizontal time
scale covers 0.000 to 0.070 s, whereas the vertical depth scale includes the
uppermost and lowermost geophone depths of 0.027 and 0.085 km, respectively.
The characteristic VSP pattern is obtained with reflections off the layer
interfaces resulting in up-going as well as down-going waves. The direct
P-wave from the source is reflected at the water-bottom interface, not only
as a P-wave but as an S-wave as well. Hence, both P- and S-waves will be
visible. They can be separated because of their different velocities giving
a less steep move-out for the S-waves. A very clear down-going S-wave is the
one that starts at the uppermost trace at about 0.007 s and reaches the
interface at z = 0.049 km at about 0.022 s.
On the uppermost trace, the Scholte interface wave is seen. It arrives at
about 0.011 s and, as expected, its amplitude is reduced drastically when
the receiver depth is increased. The strong down-going waves appearing at
times about 0.035 s on the uppermost trace are due to interfering
reflections from the ice layer and the free surface, causing a distorsion of
the pulse shape.
The matching depth was chosen at the water bottom (z = 0.025 km), but this
particular choice is not critical. However, it is critical indeed that
stabilization be performed during the back-propagation step (section 5 in
Ref. 6). Without stabilization, spurious growth would appear in the z
direction causing overflow and program breakdown. Of course, it would not
be possible to reproduce the evanescent character of the Scholte wave.
In this example, with densely spaced receivers, the modified stabilization
procedure of Appendix A in Ref. 6 is not needed. It has its main merit for
cases with thick homogeneous layers that do not contain receivers. For
example, if we were interested in the output at the lowermost geophone only,
or at a geophone at the ice boundary, this procedure would enhance
efficiency significantly.
The input to rpressfw for this VSP example was:
2 ! kfwtyp
1 ! iystab
1024 1 100000 ! lx,li(top),icworkdim
vspfig.mfw ! film
1 1.0 0 ! mhq0,freqmhq0, iazimi
2 ! icuttyp
1D-5 1D-7 .1 1.05 .0 1.0D10 ! epsglobln,epsglob,epsfakt,hfaktmax,hoff,eps..
1D-7 -0.25 ! epsfsub,hstepfsub
0.00001 2 0 2 5 ! preeps0, ixtyp0,nrus0pre,irus0,jrus0
0.0 0.070 0.0 ! tmin,dauer,ured
0.01 ! preala
1 ! iurot
500.0 ! freqdom
0.010 0.000 ! rmin,rint
vspfig.dat ! filo
vspfigtr ! maintr
! filg
0 ! intrest
1 ! ifilonpre
% ! bessfil
-1 ! iadintxfi
1.0 0 8 ! smaxpre,nelmaxpre,nbackmax
2 9 2 7 ! ixtyp,nrus,irus,jrus
0.0 20.0 430.0 ! argdiflim1,argdiflim0,argdiflim2
1 ! nparts
2 2 ! ibessarr(),npuarr()
(0.0,0.0) (1.5,0.0)
0.002 0.95 ! hmaxdh,hmaxrelwp
0.00 0.0001 ! epsh01q,preeps
0.00 0 1.0 ! hkvarpart,lunkvar,faktkvar
1 ! iconttyp
2 2 ! itail(1),itail(2)
-2 45.0 ( 1.5,0.0) ! ibesstail(1) argrktail(1) up0tail(1)
-2 -45.0 (-1.5,0.0) ! ibesstail(2) argrktail(2) up0tail(2)
7.0 ! anphpminterm
0.00001 ! preepsterm
1 0.0001 ! ircont,preepscont
The medium infile vspfig.mfw had the following appearance:
0.001 0.025 1 1 ! matching depth 0.025 km
1
3 (1.0d0,0.0d0)
0.000 0.000 0.000 0.0 0.0
0.000 0 0
3.500 1.600 0.980 0.100 0.100 0 0
0.001 0 0
3.500 1.600 0.980 0.100 0.100 1 0
0.001 0 0
1.500 0.000 1.000 0.000 0.000 0 0
0.025 0 0
1.500 0.000 1.000 0.000 0.000 1 0
0.025 0 0
2.000 1.200 2.200 0.075 0.150 0 0
0.027 0 0
2.000 1.200 2.200 0.075 0.150 2 1
0.049 0 0
2.000 1.200 2.200 0.075 0.150 11 11
0.049 0 0
4.000 2.300 2.650 0.060 0.120 0 0
0.061 0 0
4.000 2.300 2.650 0.060 0.120 6 6
0.061 0 0
3.000 1.600 2.400 0.067 0.130 0 0
0.075 0 0
3.000 1.600 2.400 0.067 0.130 7 7
0.075 0 0
6.500 4.000 3.200 0.0505 0.10 0 0
0.085 0 0
6.500 4.000 3.200 0.0505 0.10 5 5
0.085 0 0
6.500 4.000 3.200 0.0505 0.10 0 0
Note the specification of the integration path with consistent "uphoern(,)"
and "up0tail(1:2)"; cf. equation (16) and Fig. 3 in Ref. 1, although there
is a sign change and a change of direction for the second tail there. The
maximum medium slowness is 1/1.200 = 0.833 s/km, and the breakaway points
+-1.5 s/km should be sufficiently far out.
Since the parameter FREQDOM is positive, the standard sinc-type
moment-tensor wavelet is used. The initial P-arrival wavelet that is seen
in the figure of the vspfig.ps postscript file is affected to some extent by
interference from the water-bottom reflection, but at depths around 0.040 km
it looks like a sign-changed derivative of a sinc-function. This is
consistent with the displacement-response formula given in connection with
the "FILM format" explanation. At the distances and dominant frequency
considered, the second term there (a "derivative term") will dominate, and
the output is sign-changed by convention. Recall also that the positive z
direction is downwards, hence outgoing and away from the source with the
source-receiver geometry in this example. Furthermore, a deflection upwards
for a trace in the figure signifies positive output.
The space between two traces in the figure corresponds to 9.91435*(10**6)/31
= 0.319818*(10**6) km for the vertical displacement. As mentioned in the
introductory section, the results would have to be scaled down to a
realistic source strength for a quantitative interpretation.
References
----------
See the information file read.me in this directory rpress.