Choose suitable parameters for the FD run based on the model requirements and hardware capabilities. The centre frequency of the input wavelet, and maximum and minimum velocities in the input model will control the grid spacing and time step requirements.
Model size
Model size is a significant concern given the expected run times for models. Compilation of the programs should be optimised for the platform on which it is running; this has now been done (attempted?) for craton at UBC (Amor and Gorman 20.3.00) but maximum model sizes have not yet been determined. We successfully ran a model that was 540 km by 100 km with a 100 m grid length which required just over 1 Gb total arrays and about 11 hours to run.
Typical models created by Tim Henstock for Deep Probe modelling used a grid spacing of 100 m and a 1.75 Hz centre frequency (0.8-3.5 Hz bandwidth). On tienshan, a four-processor computer at Rice, these parameters enable a model to be run which is up to 1949 km long and 138 km deep. Even though tienshan (the four-processor computer at Rice) has 4 Gb memory, the compiler will only allow up to 1.75 Gb total arrays. (Tim tried from February to June to get the new compiler which will allow all the memory to be used but without any success!)
Free surface
VE2D22 cannot accommodate structure at the surface; this has been attempted in later versions of the code, but these enhancements are not part of this version (2.2). The top of the input model must be a line. For crustal modelling, it may be desirable to start the modelling at a datum beneath the surface, perhaps even beneath the sedimentary layer if one is present. From a computational viewpoint, this will have the added benefit of increasing the value of the lowest velocity in the model which will increase the maximum freqency that can be included in the forward model without dispersion. However, it should be noted that stripping off the near surface can have an adverse effect on surface wave attenuation which will result in very large surface waves relative to other phases in the model.
Rules of thumb for selecting dt and fmax (from Alan Levander)
The following two rules are meant to ensure non-dispersive propagation of the input waveform in a forth-order finite-difference routine.
dt =< 0.75*sqrt(3/8)*(h/vmax) where h is the dimension of a single cell in the grid (they must be square cells)
fmax = vmin / (5*h)
To be absolutely free of grid dispersion, the 5 in the last formula should be 10.
Radial Earth models
If a radial Earth RAYINVR model has been used, then an Earth flattening transform should be performed on the model in order to approximate a flat Earth representation of the data. A routine has been produced for this which converts a radial Earth v.in file to a flat Earth v.in file (see Andrew Gorman).
In this step, start with a v.out file from RAYINVR with the model you want to use. This v.out file must have the same sample spacing in the x-direction as it does in the z-direction. This sample spacing does not need to be as fine as that which you choose for the FD model, but it should be small enough so that little interpolation is required; it may be desireable to use a multiple of the desired FD model spacing.
This file will be modified (see next paragraph) to correspond exactly with the FD model domain to be used for input to VE2D22. It will be windowed and shifted to start at 0 depth and 0 offset. (Starting at 0 offset is not strictly necessary, but later input will require this conversion anyway, so it is convenient to do it here.) The modified model can only contain positive offset distances in the x direction. Since the code does not handle topography (see above), you need to choose a sensible datum.
Modifying v.out file
Utility: conv_vout
Extract a subset of the model using conv_vout (a modification of a routine, conv_vin from Tim Henstock). This routine will ask for four values (xmin, xmax, zmin and zmax) and will create a windowed file, v.out.new.
IMPORTANT:Two minor edits are required to v.out.new. The declared size and depth of the model (line 2 of v.out.new) must be reduced by one grid point. (Tim says that, "basically there is a bug in a later piece of code put in to cope with a bug in rayinvr....")
3. Creation of Acoustic Model
Utility: setupa3
Once the desired model exists in a v.out file, run the utility program setupa3. This requires an input control file to be piped in to it (see notes in FD Utility section.)
This produces a file named MTRPRP which is an acoustic model file.
Cell size
Cell size is a tricky thing. setupa3 will fill out the model by interpolating to the desired cell size. It will also fill 9.999's with more sensible values. It is not necessary to create a RAYINVR v.out with a small (e.g., 100 m) output spacing (which takes forever - and may be limited by space requirements of the system!). However, it is important to consider that the model is going to interpolated to create the FD input model, so the v.out file should have a sample spacing that is in the order of a km or so.
Viewing the FD acoustic model
The utility routine ac2gmt will create a three column output file of x, z and v which can then be read, gridded and countoured by gmt (or other packages) to produce a colour image of the FD model.
(At Rice, Alan has a program called rdamod to generate a segy image of the acoustic model - Tim prefers to check the visco-elastic model after the next step.)
4. Creation of Visco-elastic Model
Utility: convert_ac_vefd2
Run the utility program convert_ac_vefd2 which uses a piped input file to interpolate Vp/Vs ratios, density, Qp and
Qs for the viscoelastic code.
This produces a (large) file MTRPRP.VE. You need to rename MTRPRP (the acoustic file) and then rename MTRPRP.VE to MTRPRP.
Viewing the FD visco-elastic model
The utility routine rdpmod_segy will produce a SEG-Y file of the input which is worth checking (in ProMAX for example) to make sure there are no glitches.
5. Execute finite difference code
Make a command file (e.g., vefd.in) for the VE code as shown below.
Once set up, very few of these parameters will need to be changed. Some of the more important, model-specific, parameters are highlighted in red.
Example vefd.in file
24 lines of single quantity input. (Comments should not be included.)
vefd.in
| Example Value |
Description |
|
| 1.75 | Central frequency of wavelet. (If Ricker wavelet is chosen (see below), this in the central frequency of a 2 octave wavelet.) |
| 10000. | Width of absorbing boundary (in m) |
| 5000. | Width of taper (in m) |
| 2. | Qp in boundary |
| 2. | Qs in boundary |
| 1 | Do you want to read all material propeties from Alan Levander's MTRPRP file? (yes=1) |
| 5.0000000000E-03 | Time step (in s) |
| 0 | Output material properties on .cfl files (1=yes/0=no) ? |
| 29500 | Number of time steps. (Hint: Calculate for farthest offset for a particular shot to keep this number as small as possible.) |
| 30000 | Time steps between snap shots. (Set in connection with previous value.) |
| 4 | Plane IC (1), point IC (2), BC (3), or
insertion of point source in one point (4)? |
| 1 | Ricker wavelet (1), impulsive CW (2),
time-series sepecified on file (3), or
an integrated Ricker wavelet (4)? |
| 2 | Quantity to output in seismogram?
(0=pressure (acoustic medium), 1=P/S energy,)
(2=horizontal/vertical component of velocity) |
| .016 | Time step between samples in seismogram (in sec) |
| 1 | How many HLA? (Horizontal Line Array, i.e., a receiver array) |
| 0 | How many VLA? (Vertical Line Array, i.e., for a vertical seismic profile) |
| 0. | z co-ordinate of HLA 1* (in m) |
| 0. | x1co-ordinate of HLA 1* (in m from left side of model.) |
| 750000. | x2 co-ordinate of HLA 1* (in m from left side of model.) |
| 500. | dx for HLA 1* (in m); must be a multiple of the bin size, h. |
| 1000. | Source: Amplitude |
| 163000. | Source: x-centre (in m from left side of model.) |
| 40. | Source: z-centre (in m) |
| -.20 | Source: Time shift (should be < 0) |
* repeated for each HLA and/or VLA designated.
Utilising multi-processors
Before execution of the routine, ensure that the parallel environment is set to the appropriate level. For craton (a two-processor machine) type:
setenv PARALLEL 2
then submit the job:
ve2d22 < vefd.in
and go for coffee.
6. Output
Utility: seistrans_segy
The program will produce two files, seismogram_ver and seismogram_hor;
if you make a symbolic link
seismogram->seismogram_ver
then run seistrans_segy, the output of the vefd code is demultiplexed into a segy file (no parameters needed - output is hse01.sgy) You can run it quite happily while the vefd code is running to check on progress - worth doing after a couple of seconds are generated.
