RUN\_MULTIPLE\_SHOTS : run multiple shots defined in SOURCE\_FILE (yes=1)\\

Note that Y denotes the vertical direction!\\

The following source types are availabe: explosive sources that excite compressional waves

only (SOURCE\_TYPE=1), point forces in x-, y-, and z-direction (SOURCE\_TYPE=2,3,4), point forces in an arbitrary direction (SOURCE\_TYPE=5) and a double couple source (SOURCE\_TYPE=6).

When using the earthquake double couple source, the orientation of the slip vector on the fault plane is defined by the angles STR, DIP and RAKE. For their definition see Figure \ref{Fig:EQS}. The earthquake source is implemented via a moment-tensor description \cite{graves:96}, while the moment-tensor components are derived from the strike, dip and rake angles \cite{aki:80}. The radiated energy of the earthquake is defined by the seismic moment, which is denoted by AMON. When using the earthquake double couple source output seismograms are in units of meters per second.

\caption{The orientation of a simplified earthquake rupture plane is described by the strike angle STR, which is measured clockwise from north to the fault line, and the angle DIP, which denotes the angle between fault plane and the horizontal plane. The slip direction on the fault plane is denoted by the angle RAKE, which is measured counterclockwise from the horizontal fault line.}

\label{Fig:EQS}

\end{figure}

Three built-in wavelets of the seismic source are available. The corresponding time functions are defined in src/wavelet.c. You may modify the time functions in this file and recompile to include your

own analytical wavelet or to modify the shape of the built-in wavelets.

...

...

@@ -899,9 +870,11 @@ your source wavelet in ASCII-format in SIGNAL\_FILE. SIGNAL\_FILE should contain

The time interval between the samples must equal the time step interval (DT) of the FD simulation (see above) ! It it therefore mostly necessary to resample/interpolate a given source time function with a smaller sample rate. You may use the matlab script mfiles/resamp.m to resample your external source signal to the required sampling interval.

The following source types are availabe: explosive sources that excite compressional waves only (SOURCE\_TYPE=1), and point forces in the x-, y-, and z-direction (SOURCE\_TYPE=2,3,4).

The force sources excite both P- and S-waves. The explosive source is located at the same position as the diagonal elements of the stress tensor, i.e. at (i+1/2, j+1/2, k+1/2) (Figure \ref{fig_cell}).

The forces are located at the same position as the corresponding components of particle velocity (Figure \ref{fig_cell}). If (x,y,z) denotes the position at which the source location is defined in source.dat, then the actual force in x-direction is located at (x+DX/2,y,z), the actual force in y-direction is located at (x,y+DY/2,z) and the actual force in z-direction is at (x,y,z+DZ/2).

The following source types are availabe: explosive sources that excite compressional waves

only (SOURCE\_TYPE=1), point forces in x-, y-, and z-direction (SOURCE\_TYPE=2,3,4), point forces in an arbitrary direction (SOURCE\_TYPE=5) and a earthquake double couple source (SOURCE\_TYPE=6).

The explosive source (SOURCE\_TYPE=1) is located at the same position as the diagonal elements of the stress tensor, i.e. at (i, j, k) (Figure \ref{fig_cell}).

The force sources (SOURCE\_TYPE=2,3,4) are located at the same position as the corresponding components of particle velocity (Figure \ref{fig_cell}). If (x,y,z) denotes the position at which the source location is defined in source.dat, then the actual force in x-direction is located at (x+DX/2,y,z), the actual force in y-direction is located at (x,y+DY/2,z) and the actual force in z-direction is at (x,y,z+DZ/2).

By choosing SOURCE\_TYPE=5 a custom directive force can be defined by the force angle between horizontal x- and z-coordinate (SOURCE\_ALPHA) and the force angle between horizontal x- and vertical y-coordinate (SOURCE\_BETA). The orientation of the custom force source angles is plotted in Figure \ref{coordinate_system_customforcesource}. The unit of both SOURCE\_ALPHA and SOURCE\_BETA is in degree from 0 to 360. Please note, that due to numerical imprecision of the calculation of the values of the sine or cosine function, there is a small difference in the seismogram output when considering, e.g., the vertical y-force source (SOURCE\_TYPE=3) and the equivalent custom force source (SOURCE\_TYPE=5, SOURCE\_ALPHA=0, SOURCE\_BETA=90). Usually, this difference is orders of magnitude smaller than the actual seismogram.

\begin{figure}

...

...

@@ -913,6 +886,38 @@ By choosing SOURCE\_TYPE=5 a custom directive force can be defined by the force

\end{figure}

When using the earthquake double couple source (SOURCE\_TYPE=6), the orientation of the slip vector on the fault plane is defined by the angles STR, DIP and RAKE. For their definition see Figure \ref{Fig:EQS}. The earthquake source is implemented via a moment-tensor description \cite{graves:96}, while the moment-tensor components are derived from the strike, dip and rake angles \cite{aki:80}. The radiated energy of the earthquake is defined by the seismic moment, which is denoted by AMON. When using the earthquake double couple source output seismograms are in units of meters per second.

\caption{The orientation of a simplified earthquake rupture plane is described by the strike angle STR, which is measured clockwise from north to the fault line, and the angle DIP, which denotes the angle between fault plane and the horizontal plane. The slip direction on the fault plane is denoted by the angle RAKE, which is measured counterclockwise from the horizontal fault line.}

\label{Fig:EQS}

\end{figure}

\clearpage

The locations of multiple sources must be defined in an external ASCII file (SOURCE\_FILE) that has the following format:

\begin{verbatim}

% XSRC YSRC ZSRC TD FC AMP

...

...

@@ -922,12 +927,12 @@ More than one source position may be specified for example to simulate acoustic

XSRC is the x-coordinate of a source point (in meter), YSRC is the vertical y-coordinate of a source point (in meter), ZSRC is the horizontal z-coordinate of a source point (in meter).

TD is the excitation time (time-delay) for the source point [in seconds], FC is the center frequency of the source signal [in Hz], and AMP is the maximum amplitude of the source signal.

Please note that SOFI3D considers real physical units as input and output. AMP specifies a seismic moment of 1 Nm if an explosive source is chosen, in case of force source AMP specifies a force of 1 N. The seismogram output scales to this seismic moment or force and is given in particle velocities or pressure. Only the component \textit{div} and \textit{curl} cannot be treated as physical measures. For more information on the seismogram output see section \subsubsection{Receivers}.

Please note that SOFI3D considers real physical units as input and output. AMP specifies a seismic moment of 1 Nm if an explosive source is chosen, in case of force source AMP specifies a force of 1 N. The seismogram output scales to this seismic moment or force and is given in particle velocities or pressure. Only the component \textit{div} and \textit{curl} cannot be treated as physical measures. For more information on the seismogram output see section \ref{Receivers}.

The SOURCE\_FILE = \lstinline{./sources/source.dat} that defines an explosive source at $x_s=2610.0\;$ m, $y_s=2610.0\;$ m and $z_s=2100.0\;$ m (depth) with

a center frequency of 5 Hz (no time delay) reads, for example

The SOURCE\_FILE = \lstinline{./sources/source.dat} that defines an explosive source at $x_s=2592.0\;$ m, $y_s=2106.0\;$ m (depth) and $z_s=2592.0\;$ m with no time delay,

a center frequency of 5 Hz and an maximum amplitude of 1.0 reads, for example

\begin{verbatim}

2592.0 2106.0 2106.0 0.0 5.0 1.0

2592.0 2106.0 2592.0 0.0 5.0 1.0

\end{verbatim}

If the option RUN\_MULTIPLE\_SHOTS=0 in the parameter file all the shot points defined in the SOURCE\_FILE are used in one simulation. Setting RUN\_MULTIPLE\_SHOTS=1 will start model runs

...

...

@@ -953,7 +958,7 @@ If READMOD=1, the P-wave, S-wave, and density model grids are read from external

In these files, each material parameter value must be saved as 32 bit (4 byte) native float. Velocities must be in meter/second, density values in $kg/m^3$. The fast dimension is the y direction. See src/readmod.c. The number of samples for the entire model in the x-direction is NX, the number of values in the y-direction is always NY and the number of values in the z-direction is always NZ. The file size of each model file thus must be NX*NY*NZ*4 bytes. You may check the model structure using the SU command ximage:

ximage n1=$<$NX$>$ n2=$<$NY$>$$<$ model/test.vp .

ximage n1=$<$NY$>$ n2=$<$NX$>$$<$ model/test.vp .

If the variable TAU = 0.0 (see the next section \textit{Q-approximation}), it is also possible to read Qp, and Qs grid files to allow for spatial variable attenuation. In case of TAU > 0, constant damping is assumed with $Q_p=Q_s =\frac{2}{\mbox{TAU}}$. Please not that in order to ensure proper visualization using xmovie or ximage, the binary models are structured in terms of coordinate priority as follows : [Y,X,Z]. That means a 3-D model consits of NZ planes of the size NY $\times$ NX. The fast dimension of each NY-NX plane is the vertical Y-direction, i.e., first, columns of NY material parameters are written to file, in total NX columns of the length NY are written. To give an example of how to create a 3-D model using Matlab, the M-File \lstinline{create_simple_tunnel_mod3D.m} is located in the folder \lstinline{mfiles} (see Chapter \ref{installation}).

...

...

@@ -1055,6 +1060,7 @@ SNAP\_PLANE : output of snapshots as energy (without sign=1, with sign true for

If SNAP$>0$, wavefield information (particle velocities, pressure, or curl and divergence of particle velocities) for the entire model is saved on the hard disk (assure that enough free space is on disk!). Each PE is writing his sub-volume to disk. The filenames have the basic filename SNAP\_FILE plus an extension that indicates the PE number in the logical processor array (see Figure \ref{fig_grid}), i.e. the PE with number PEno writes his wavefield to SNAPFILE.PEno. The first snapshot is written at TSNAP1 seconds of seismic wave traveltime to the output files, the second at TSNAP1+TSNAPINC seconds etc. The last snapshots contains wavefield at TSNAP2 seconds. Note that the file sizes increase during the simulation. The snapshot files might become quite LARGE. It may therefore be necessary to reduce the amount of snapshot data by increasing IDX, IDY and IDZ and/or TSNAPINC. A detailed description how to visualize 3-D wavefields is given in section \ref{visual}. In order to merge the separate snapshot of each PE after the comletion of the wave modeling, you can use the program snapmerge (see Chapter \ref{installation}, section \textbf{src}). The bash command line to merge the snapshot files can look like this: \lstinline{../bin/snapmerge ./in_and_out/sofi3D.json}.

\subsubsection{Receivers}

\label{Receivers}

\begin{verbatim}

"Receiver" : "comment",

"SEISMO" : "4",

...

...

@@ -1079,7 +1085,7 @@ NGEOPH : distance between two adjacent receivers (in gridpoints)\\

If SEISMO$>$0, seismograms are saved on hard disk. If SEISMO equals 1 x-, y-, and z-component of particle velocity will be written according to parameters specified in Chapter \ref{seismograms}.

If SEISMO==2 pressure (sum of the diagonal components of the stress tensor) recorded at the receiver locations (receivers are hydrophones!) is written. if SEISMO=3 the curl and divergence are saved. As described in section \subsubsection{Sources} the seismogram output is physically correct for a seismic moment of 1 Nm (explosive source) and a force of 1 N (force source), respectively. A validation of SOFI3D and analytical data can be found in section \ref{analytic}.

If SEISMO==2 pressure (sum of the diagonal components of the stress tensor) recorded at the receiver locations (receivers are hydrophones!) is written. if SEISMO=3 the curl and divergence are saved. As described in section \ref{Sources} the seismogram output is physically correct for a seismic moment of 1 Nm (explosive source) and a force of 1 N (force source), respectively. A validation of SOFI3D and analytical data can be found in section \ref{analytic}.

The curl and divergence of the particle velocities are useful to separate between P- and S-waves in the snapshots of the wavefield. SOFI3D calculates the divergence and the magnitude of the curl of the particle velocity field \cite{dougherty:88}. The motivation for this is as follows. According to Morse and Feshbach \shortcite{morse:53} the energy of P- and S-wave particle velocities is, respectively,

\begin{equation}

...

...

@@ -1331,7 +1337,7 @@ The command to start a simulation on 8 processor with the lowest priority of -19

If you use LAMMPI the command \lstinline{lamboot -v lamhost} must be executed on node 0 which is the PE where ./par/lamhosts is the file containing IP addresses of all computing nodes. It if often useful to save the standard output of the program for later reference. The screen output may be saved to sofi3D.jout using

Combining both the screen and the file output you can execute the following line (this is used in the shell script \lstinline{startSOFI3D.sh})

...

...

@@ -1506,7 +1512,8 @@ reflections at the boundaries of the model grid.

In order to gain an absolute validation of Finite Difference data gained from SOFI3D, we perform comparisons of SOFI3D seimogram traces and data calculated by analytic solutions. The analytical solution of a point explosive source within a homogeneous full space and a point vertical force source at the top of a homogenous halfspace are provided by the Matlab Script \lstinline{analytic_solution_main.m} located in the folder \lstinline{mfiles/analytical_solution_fullspace} and \lstinline{mfiles/analytical_solution_halfspace}. In each folder you also find the Matlab script \lstinline{comparison_analytic_vs_sofi3D.m} to directly compare the SOFI3D and analytical data and gain the same figures which are presented in the following.

We first consider a point explosive source within a homogeneous full space. We find this validation so useful that we already included the SOFI3D modeling into the overnight built test. For further information on the input parameters and elastic properties we refer to the Documentation of the SOFI3D overnight built (see \lstinline{doc/overnightbuilt_sofi3D/sofi3D_overnight_documentation.pdf}, section \textit{Homogeneous Fullspace}). For several offsets the SOFI3D and analytical traces are displayed in Figure \ref{analytic_sofi3D_p} (pressure component) and Figure \ref{analytic_sofi3D_vx} ($v_x$ component). Each Figure contains the analytical data (red line), SOFI3D data (dashed blue line) and the difference (dotted green line). Both the FD and the analytical data are tracewise normalized to the maximum of the analytical data. The amplitudes are thus truly scaled and the difference is absolute. For each trace the RMS (L2 Norm) is calculated and the RMS error given by

We first consider a point explosive source within a homogeneous full space. We find this validation so useful that we already included the SOFI3D modeling into the overnight built test. For further information on the input parameters and elastic properties we refer to the Documentation of the SOFI3D overnight built

(see \lstinline{doc/overnightbuilt_sofi3D/sofi3D_overnight _documentation.pdf}, section \textit{Homogeneous Fullspace}). For several offsets the SOFI3D and analytical traces are displayed in Figure \ref{analytic_sofi3D_p} (pressure component) and Figure \ref{analytic_sofi3D_vx} ($v_x$ component). Each Figure contains the analytical data (red line), SOFI3D data (dashed blue line) and the difference (dotted green line). Both the FD and the analytical data are tracewise normalized to the maximum of the analytical data. The amplitudes are thus truly scaled and the difference is absolute. For each trace the RMS (L2 Norm) is calculated and the RMS error given by