**# Introduction**
Elastic Lattice Method (ELM) is a discrete particle numerical method for simulating the seismic waves propagation through complex media in the presence of topography. The structural model consists of particles arranged on a cubic lattice which interact through a central force term and bond-bending force. Particle disturbances are followed through space by numerically solving their equation of motion. The description of the code, examples, and benchmarking against finite difference codes can be found in [1-6] . The ELM code has both 2D and 3D implementations, but the input files are nearly identical. Here we take the 2D version for instruction.

**# Potential application in volcano observatories**
The aim of volcano monitoring is to understand physical processes within a volcano and ultimately help in forecasting potential hazards. One of the primary interests in volcano-seismology is determining volcanic seismic sources, but the generation and propagation of seismic waves throughout volcanic regions is a complex and nonlinear phenomenon controlled by the interaction of many processes. Because volcanoes are highly mechanically heterogeneous they significantly distort seismic wave propagation and hence the source ‘fingerprint’ in seismograms is often hidden by ‘path effects’ created when the wave propagates from to the source to the seismic receiver. Using the VolcSeisSimu tool researchers or monitoring scientists can determine these propagation path effects for arbitrarily heterogeneous models and ‘deconvolve’ them from observed seismograms to better understand and constrain volcano source models. Source inversion is a highly specialised and time consuming research area, but is critical to a better understanding of volcano seismograms in terms of source processes. This scheme offers an alternative 3D method for modelling wave propagation in the presence of strong topography and subsurface heterogeneity, and the flexibility of the method allows for a wide range of possible source mechanisms with no restriction on source geometric shape or on the distribution and degree of elastic heterogeneity.

**# Getting Started**
To download the ELM, both of 2D and 3D, software packages, type this:

`git clone --recursive https://git.dias.ie/graphiit/DIASGEO.git`

Or if you prefer, download the source code directly in zip or tar format by clicking on the download button on this page:

https://git.dias.ie/graphiit/DIASGEO

The program is written in C with MPI commands for parallel computing. It consists of the main code `ELM.c`

, input parameter file `PARAMETERS.h`

, and a number of subroutines. All these files are located within a single directory `src`

. Please keep in mind that when you change the head file `PARAMETERS.h`

, you should redo the compilation before running.

To compile the ELM, type this:

`mpicc ELM.c -o ELM -lm`

To execute the ELM, type this:

`mpiexec -np nproc ELM`

where `nproc`

is the number of processes.

You can also simply run the Bash script `run_example.sh`

to compile and execute ELM as the one-button-solution.

# Setting the inputs

An example of `PARAMETERS.h`

file is shown below, followed by detailed explanations of every line/section of the file.

```
//---------------------------------------------------------------------------
// PARAMETERS.h
//---------------------------------------------------------------------------
// MODEL VALUE
#define Max_Time (8001) // No of time steps
#define Nx ( 55) // Length of x grid (mpi)
#define Ny ( 495 ) // Length of y grid (no mpi)
//---------------------------------------------------------------------------
// BOUNDARY CONDITIONS
#define absorb_model 2 // 0 -> off 1-> full on: 2 -> surface on;
#define lambda (2.5e-6) // decay constant for boundary (1/3n^2)
#define width 100 // no of absorbing columns
#define flat_surface 0 // 0- no flat surface 1-flat surface at
#define topography 1 // 0-topography off ; 1-> topography on
#define input5 "./input/topo" // topography file (x,y)
#define velocity_model 0 // 0-homogeneous: 1-input file:
#define input4 "./input/vp" // P-velocity file
//---------------------------------------------------------------------------
// PHYSICAL VALUES
#define dt (0.0005) // Time step
#define dx (10.0) // Grid spacing
#define density (2000.0) // Density in kg/m^3
#define velp (3500.00) // P-wave in m/s
#define vels (2500.00) // S-wave in m/s
//---------------------------------------------------------------------------
// SOURCE
#define input_source 1 // 1-external input source; 3-no source 4-MT
#define input1 "./input/source" // Source location file
#define NumSource 1 // Number of entries in source_locations file
#define sdir 3 // Source direction 0=x 1=z 3=explosive
#define sscale (1e+8) // Scale source amplitude by this factor (Force)
#define input2 "./input/signal" // Source file
#define mxx 1.0
#define mzz 1.0
#define mxz 0.0
//---------------------------------------------------------------------------
// RECEIVERS && SNAPSHOTS
#define NumGeo 2 // Number of entries in geophone file
#define input0 "./input/stations" // Trace coordinates in x y format
#define NumSnap 40 // No of snapshots
#define isnap 200 // Snapshot record interval
#define snap_start 200 // Start snapshots
//---------------------------------------------------------------------------
// Fracture models
#define fracture_model 0 // Fractures on 1 or off 0
#define NoFrac 105 // No of fractures
#define input3 "fractures" // Fracture file
//---------------------------------------------------------------------------
// I/O FILES
#define output2 "vel" // Velocity snapshots
#define output6 "disp" // Displacement snapshots
#define output4 "station" // Seismograms
//---------------------------------------------------------------------------
```

The specification of `MODEL VALUE`

:

`Max_time`

: the total number of time steps.

`Nx`

: Number of data points in x direction, which will be run on each node. In this example, if the simulation is run on `nproc`

nodes so the total number of grid points in x direction is `nproc x Nx`

.

`Ny`

: number of grid points in y direction.

The specification of `BOUNDARY CONDITIONS`

:

In order to make the boundaries of the numerical domain “transparent” for the seismic waves (i.e. without artificial reflections from the edges of the model), absorbing boundaries are put around the model. The equation for the absorbing boundary is of the form:

ω=e^{-λx},

that the wavefield is slowly damped during the propagation through the absorbing boundaries and back. The choice of the thickness of the boundaries x (in the number of grid points) and decay constant λ are of crucial importance to suppress unwanted reflection.

`absorb_model`

: `0`

switch off all absorbing boundaries; `1`

switch on all absorbing boundaries; `2`

all absorbing boundaries except at the free surface are switched on.

`lambda`

: decay constant.

`width`

: thickness of absorbing boundaries in grid points

`topography`

: `0`

don't use topography, that is flat top; `1`

use topography.

`input5`

: specify topography file. Each line contains 2 values (x, topo), where x=0, 1, .... Ensure that the topography file does not have any single points sticking out above their neighbouring points, it can cause high-frequency unrealistic resonance of your numerical model.

`velocity_model`

: `0`

homogeneous model, and `1`

heterogeneous model.

`input4`

: heterogeneous model file with the P-wave velocity specification; each line contains 4 values (x, y, vp), where x=0, 1, ..., Nx-1, and y=0, 1, ..., Ny-1. Currently, the S-wave velocity and density are calculated from the P-wave velocity inside the code with fixed Poisson’s ratio of 1/4, that

vs = vp/sqrt(3) and ρ= 1700+0.2vp,

but the value of Possion's ratio can be changed as desired.

The specification of `PHYSICAL VALUES`

:

`dt`

: time step used in the simulation.

From the stability criterion, the time step should satisfies

dt<dx/(2v_{max}),

where v_{max} is the maximal velocity of the model.

`dx`

: grid spacing in meters (the grid spacing is equal in all directions).

From the dispersion criterion, there should be more than 25 spatial points to sample the minimal wavelength, i.e.,

dx < λ_{min}/25,

for the model including topography.

`density`

: density value in homogeneous model.

`velp`

: p-wave velocity value in homogeneous model.

`vels`

: s-wave velocity value in homogeneous model.

The specification of `SOURCE`

:

`input_source`

: `0`

without source; `1`

with source.

`NumSource`

: specify the number of simultaneous sources you would like to use.

`sdir`

: `0`

horizontal single force in x direction; `1`

single force in y direction; `3`

a purely isotropic source; `4`

a moment tensor source.

`sscale`

: factor for scaling the source-time function.

`input3`

: the source wavelet file name, which is a single column file.

`input2`

: the source location file name.

`mxx`

, `mzz`

, `mxz`

: define the moment tensor matrix.

The specification of `RECEIVERS AND SNAPSHOTS`

.

`NumGeo`

: number of receivers.

`input0`

: name of the file containing location coordinates for all receivers.

`NumSnap`

: number of volume snapshots.

`isnap`

: number of time steps between two consecutive snapshots.

`snap_start`

: number of time steps of the first snapshot.

The specification of `OUTPUT FILE NAMES`

:

`output2`

: name prefix for velocity snapshots.

`output4`

: name prefix for seismograms.

`output6`

: name prefix for displacement snapshots.

In summary, in addition to the `PARAMETERS.h`

file, there are generally another 5 files as the input to be set:

- station locations
- source wavelet
- source location
- velocity structure
- topography

And their file names are specified in `PARAMETERS.h`

.

# Examples

Two ELM running examples with the associated pre-processing Octave/Matlab scripts and post-processing GMT scripts are provided in the folder named `elm2d`

. The execution of ELM is scripted with the Bash code `run_example.sh`

.

Example 1:

In the first example, a model with Gaussian top topography and gradient velocity is built, and the source is selected as the explosive source in the middle. The seismic wavefield snapshots are taken and the seismograms are recorded by two stations symmetrically located on the free surface.

Fig 1-1: The model deployment.

Fig 1-2: Animation of seismic wavefield.

Fig 1-3: Seismograms.

Example 2:

In the second example, 3 different type of sources are taken in the homogeneous model, e.g., a single force in X direction, an explosive source, and a dip-slip moment tensor source. The seismic wavefield snapshots are taken at the same time.

Fig 2: Snapshots of seismic wavefield.

# References

[1] O'Brien, Gareth S., Bean, Christopher J. (2004), A 3D discrete numerical elastic lattice method for seismic wave propagation in heterogeneous media with topography, Geophys. Res. Lett., Vol. 31, No. 14, L14608 10.1029/2004GL020069.

[2] O'Brien, G. S., Bean, C. J. (2004), A discrete numerical method for modeling volcanic earthquake source mechanisms, J. Geophys. Res., Vol. 109, No. B9, B09301 10.1029/2004JB003023.

[3] O'Brien, G. S. (2008), Discrete visco-elastic lattice methods for seismic wave propagation, Geophys. Res. Lett., 35, L02302, doi:10.1029/2007GL032214.

[4] O'Brien G.S., Bean C.J. and Tapamo H. (2009), Dispersion analysis and computational efficiency of Elastic Lattice methods for seismic wave propagation, Computers & Geosciences, Vol. 35, 1768-1775.

[5] O'Brien, G. S. and C.J. Bean (2011), An irregular lattice method for elastic wave propagation, Geophysical Journal International, doi: 10.1111/j.1365-246X.2011.05229.x.

[6] O'Brien, G.S., T. Nissen-Meyer and C.J. Bean (2012), A Lattice Boltzmann method for elastic wave propogation in a Poisson solid, Bulletin of the Seismological Soceity of America, VOLUME 103, No. 3, Pages 1224-1234, doi: 10.1785/0120110191.

# Links

This software was made available as part of the EUROVOLC project. For more information see https://eurovolc.cp.dias.ie/index.php/Open_software