09_noise_simulations.tex

\chapter{Noise Cross-correlation Simulations}


Besides earthquake simulations, SPECFEM3D Cartesian includes functionality
for seismic noise tomography as well. In order to proceed successfully
in this chapter, it is critical that you have already familiarized
yourself with procedures for meshing (Chapter \ref{cha:Mesh-Generation}),
creating distributed databases (Chapter \ref{cha:Creating-Distributed-Databases}),
running earthquake simulations (Chapters \ref{cha:Running-the-Solver})
and adjoint simulations (Chapter \ref{cha:Adjoint-Simulations}).
Also, make sure you read the article `Noise cross-correlation sensitivity
kernels' \citep{trompetal2010}, in order to understand noise simulations
from a theoretical perspective.


\section{Input Parameter Files}

As usual, the three main input files are crucial: \texttt{Par\_file},
\texttt{CMTSOLUTION} and \texttt{STATIONS}. Unless otherwise specified,
those input files should be located in directory \texttt{DATA/}.\\


\texttt{CMTSOLUTION} is required for all simulations. At a first glance,
it may seem unexpected to have it here, since the noise simulations
should have nothing to do with the earthquake -- \texttt{CMTSOLUTION}.
However, for noise simulations, it is critical to have no earthquakes.
In other words, the moment tensor specified in \texttt{CMTSOLUTION}
must be set to zero manually!\\


\texttt{STATIONS} remains the same as in previous earthquake simulations,
except that the order of receivers listed in \texttt{STATIONS} is
now important. The order will be used to determine the `main' receiver,
i.e., the one that simultaneously cross correlates with the others.\\


\texttt{Par\_file} also requires careful attention. A parameter called
\texttt{NOISE\_TOMOGRAPHY} has been added which specifies the type
of simulation to be run. \texttt{NOISE\_TOMOGRAPHY} is an integer
with possible values 0, 1, 2 and 3. For example, when \texttt{NOISE\_TOMOGRAPHY}
equals 0, a regular earthquake simulation will be run. When it is
1/2/3, you are about to run step 1/2/3 of the noise simulations respectively.
Should you be confused by the three steps, refer to \citet{trompetal2010}
for details.\\


Another change to \texttt{Par\_file} involves the parameter \texttt{NSTEP}.
While for regular earthquake simulations this parameter specifies
the length of synthetic seismograms generated, for noise simulations
it specifies the length of the seismograms used to compute cross correlations.
The actual cross correlations are thus twice this length, i.e., $2~\mathrm{NSTEP}-1$.
The code automatically makes the modification accordingly, if \texttt{NOISE\_TOMOGRAPHY}
is not zero.\\


There are other parameters in \texttt{Par\_file} which should be given
specific values. For instance, since the first two steps for calculating
noise sensitivity kernels correspond to forward simulations, \texttt{SIMULATION\_TYPE}
must be 1 when \texttt{NOISE\_TOMOGRAPHY} equals 1 or 2. Also, we
have to reconstruct the ensemble forward wavefields in adjoint simulations,
therefore we need to set \texttt{SAVE\_FORWARD} to \texttt{.true.}
for the second step, i.e., when \texttt{NOISE\_TOMOGRAPHY} equals
2. The third step is for kernel constructions. Hence \texttt{SIMULATION\_TYPE}
should be 3, whereas \texttt{SAVE\_FORWARD} must be \texttt{.false.}.\\


Finally, for most system architectures, please make sure that \texttt{LOCAL\_PATH}
in \texttt{Par\_file} is in fact local, not globally shared. Because
we have to save the wavefields at the earth's surface at every time
step, it is quite problematic to have a globally shared \texttt{LOCAL\_PATH},
in terms of both disk storage and I/O speed.


\section{Noise Simulations: Step by Step}

Proper parameters in those parameter files are not enough for noise
simulations to run. We have more parameters to specify: for example,
the ensemble-averaged noise spectrum, the noise distribution etc.
However, since there are a few `new' files, it is better to introduce
them sequentially. In this section, standard procedures for noise
simulations are described.


\subsection{Pre-simulation}
\begin{itemize}
\item As usual, we first configure the software package using:
\begin{verbatim}
./configure FC=ifort MPIFC=mpif90
\end{verbatim}

Use the following if SCOTCH is needed:
\begin{verbatim}
./configure FC=ifort MPIFC=mpif90 --with-scotch-dir=/opt/scotch
\end{verbatim}

\item Next, we need to compile the source code using:
\begin{verbatim}
make xgenerate_databases
make xspecfem3D
\end{verbatim}

\item Before we can run noise simulations, we have to specify the noise
statistics, e.g., the ensemble-averaged noise spectrum. Matlab scripts
are provided to help you to generate the necessary file:
\begin{verbatim}
EXAMPLES/noise_tomography/NOISE_TOMOGRAPHY.m  (main program)
EXAMPLES/noise_tomography/PetersonNoiseModel.m
\end{verbatim}

In Matlab, simply run:
\begin{verbatim}
NOISE_TOMOGRAPHY(NSTEP, DT, Tmin, Tmax, NOISE_MODEL)
\end{verbatim}

\texttt{DT} is given in \texttt{Par\_file}, but {\color{red} \texttt{NSTEP} is
NOT the one specified in \texttt{Par\_file}}. {\bf Instead, you have to
feed $2~\mathrm{NSTEP}-1$ to account for the doubled length of cross
correlations}. \texttt{Tmin} and \texttt{Tmax} correspond to the period
range you are interested in, whereas \texttt{NOISE\_MODEL} denotes
the noise model you will be using
(\texttt{'NLNM'} for New Low Noise Model or \texttt{'NHNM'} for New High Noise Model).
Details can be found in the Matlab script.


After running the Matlab script, you will be given the following information
(plus a figure in Matlab):
\begin{verbatim}
*************************************************************
the source time function has been saved in:
..../S_squared (note this path must be different)
S_squared should be put into directory:
./NOISE_TOMOGRAPHY/ in the SPECFEM3D Cartesian package
\end{verbatim}

In other words, the Matlab script creates a file called \texttt{S\_squared},
which is the first `new' input file we encounter for noise simulations.

One may choose a flat noise spectrum rather than Peterson's noise
model. This can be done easily by modifying the Matlab script a little.

\item Create a new directory in the SPECFEM3D Cartesian package, name it
as \texttt{./NOISE\_TOMOGRAPHY/}. We will add some parameter
files later in this folder.

\item Put the Matlab-generated-file \texttt{S\_squared} in \texttt{./NOISE\_TOMOGRAPHY/}.

That's to say, you will have a file \texttt{./NOISE\_TOMOGRAPHY/S\_squared}
in the example provided in the SPECFEM3D Cartesian package.

\item Create a file called \texttt{./NOISE\_TOMOGRAPHY/irec\_main\_noise}.
Note that this file is located in directory \texttt{./NOISE\_TOMOGRAPHY/}
as well. In general, all noise simulation related parameter files
go into that directory. \texttt{irec\_main\_noise} contains only
one integer, which is the ID of the `main' receiver. For example,
if this file contains 5, it means that the fifth receiver listed in
\texttt{DATA/STATIONS} becomes the `main'. That's why we mentioned
previously that the order of receivers in \texttt{DATA/STATIONS} is
important.


Note that in some simulations, the \texttt{DATA/STATIONS} might
contain receivers which are outside of our computational domains.
Therefore, the integer in \texttt{irec\_main\_noise} is actually
the ID in \texttt{DATA/STATIONS\_FILTERED} (which is generated by
\texttt{bin/xgenerate\_databases}).

\item Create a file called \texttt{./NOISE\_TOMOGRAPHY/nu\_main}.
This file holds three numbers, forming a (unit) vector. It describes
which component we are cross-correlating at the `main' receiver,
i.e., ${\hat{{\bf \nu}}}^{\alpha}$ in \citet{trompetal2010}. The
three numbers correspond to E/N/Z components respectively. Most often,
the vertical component is used, and in those cases the three numbers
should be 0, 0 and 1.

\item Describe the noise direction and distributions in \texttt{src/specfem3d/noise\_tomography.f90}.
Search for a subroutine called \texttt{noise\_distribution\_direction}
in \texttt{noise\_tomography.f90}. It is actually located at the very
beginning of \texttt{noise\_tomography.f90}. The default assumes vertical
noise and a uniform distribution across the whole free surface of
the model. It should be quite self-explanatory for modifications.
Should you modify this part, you have to re-compile the source code.
\end{itemize}

\subsection{Simulations}

With all of the above done, we can finally launch our simulations.
Again, please make sure that the \texttt{LOCAL\_PATH} in \texttt{DATA/Par\_file}
is not globally shared. It is quite problematic to have a globally
shared \texttt{LOCAL\_PATH}, in terms of both disk storage and speed
of I/O (we have to save the wavefields at the earth's surface at every
time step).

As discussed in \citet{trompetal2010}, it takes three steps/simulations
to obtain one contribution of the ensemble sensitivity kernels:
\begin{itemize}
\item Step 1: simulation for generating wavefields
\begin{verbatim}
SIMULATION_TYPE = 1
NOISE_TOMOGRAPHY = 1
SAVE_FORWARD (not used, can be either .true. or .false.)
\end{verbatim}

\item Step 2: simulation for ensemble forward wavefields
\begin{verbatim}
SIMULATION_TYPE = 1
NOISE_TOMOGRAPHY = 2
SAVE_FORWARD = .true.
\end{verbatim}

\item Step 3: simulation for ensemble adjoint wavefields and sensitivity
kernels
\begin{verbatim}
SIMULATION_TYPE = 3
NOISE_TOMOGRAPHY = 3
SAVE_FORWARD = .false.
\end{verbatim}

Note Step 3 is an adjoint simulation, please refer to previous chapters
on how to prepare adjoint sources and other necessary files, as well
as how adjoint simulations work.

\end{itemize}

It is better to run the three steps continuously within
the same job on a cluster, otherwise you have to collect the saved surface movies
from the old nodes to the new nodes. This process varies
from cluster to cluster and thus cannot be discussed here. Please
ask your cluster administrator for information/configuration of the
cluster you are using.\\



\subsection{Post-simulation}

After those simulations, you have all stuff you need, either in the
\texttt{OUTPUT\_FILES/} or in the directory specified by \texttt{LOCAL\_PATH}
in \texttt{DATA/Par\_file} (which are most probably on local nodes).
Collect whatever you want from the local nodes to your workstation,
and then visualize them. This process is the same as what you may
have done for regular earthquake simulations. Refer to other chapters
if you have problems.\\


Simply speaking, two outputs are the most interesting: the simulated
ensemble cross correlations and one contribution of the ensemble sensitivity
kernels.\\


The simulated ensemble cross correlations are obtained after the second
simulation (Step 2). Seismograms in \texttt{OUTPUT\_FILES/} are actually
the simulated ensemble cross correlations. Collect them immediately
after Step 2, or the Step 3 will overwrite them. Note that we have
a `main' receiver specified by \texttt{irec\_main\_noise}, the
seismogram at one station corresponds to the cross correlation between
that station and the `main'. Since the seismograms have three components,
we may obtain cross correlations for different components as well,
not necessarily the cross correlations between vertical components.\\


One contribution of the ensemble cross-correlation sensitivity kernels
are obtained after Step 3, stored in the \texttt{DATA/LOCAL\_PATH}
on local nodes. The ensemble kernel files are named the same as regular
earthquake kernels.\\


You need to run another three simulations to get the other contribution
of the ensemble kernels, using different forward and adjoint sources
given in \citet{trompetal2010}.


\section{Example}

In order to illustrate noise simulations in an easy way, one example
is provided in \texttt{EXAMPLES/noise\_tomography/}. See \texttt{EXAMPLES/noise\_tomography/README}
for explanations. \\


Note, however, that they are created for a specific workstation (CLOVER@PRINCETON),
which has at least 4 cores with `mpif90' working properly. \\


If your workstation is suitable, you can run the example in \texttt{EXAMPLES/noise\_tomography/}
using:\\

\texttt{./pre-processing.sh}\\


Even if this script does not work on your workstation, the procedure
it describes is universal. You may review the whole process described
in the last section by following the commands in \texttt{pre-processing.sh},
which should contain enough explanations for all the commands.