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.