Sections of the GTdef model
(Clicking on section heads will return to the Table of Contents)

Manual developed by Derrick Murekezi and Andrew Newman
GTdef is an implementation of the Okada (1985) dislocation model written by
Lujia Feng, Ting Chen, and Andrew Newman, with codes translated from disl (by Peter Cervelli). Okada (1985) can describe the
surface deformation associated with rectangular dislocations of most orientations, within
a purelyelastic halfspace, however we've further adapted the codes to allow for horizontally layered materials.
Fault planes do require one edge to be horizontal. Dislocations can be described with dipslip, strikeslip, and/or opening,
and hence can be used to describe a range of processes including earthquake faulting, interseismic
coupling, and planar fluid changes like magmatic dike inflation or groundwater withdrawal.
The code is developed so that one could perform both forward models that predict ground deformation
and inverse models that use data with errors and some modeling constraints to predict information about the
source. The source in either case can be either simple rectangular slips along the plane, orthogonal to it
or some combination of both.
Requirements:
 Algorithms are written in Matlab and require the 'Parameter estimation' toolbox for inversions.
 If you wish to use parallel processors, your version needs to include parpool, which was implemented around 2015.
 For plotting GMT version 5 or 6.

Set up Instructions:
 In the selected directory, to clone GTdef run
#this will create a copy of the code which you can pull to every time there is an update to the master git repository
git clone https://github.com/avnewman/GTDef.git
 Be sure the current version of the programs are in your path. You can check this using the path function, and add it using addpath.
# within matlab
> addpath('/path/to/GTdef/v4.0/stable/') %replacing with your actual path
 Additiionally, create a separate directory to run examples and save files. Running Matlab calls in this directory will save your results in one place to be analysed later.

Kinematic Modeling:
A very brief discussion of kinematic modeling methods used here should exist. This would include each Okada (1985) analytical models,
the Wang et al (2003) layered method, and general injestion of Green's functions from other sources.
Inverse Modeling:
GTdef performs a bounded linear leastsquares inversion of Green's functions to define
the bestfitting set of models to describe a dataset with errors. The program uses the parameter optimization toolbox program called lsqlin (Matlab info on it).
You can read more about it and general inversion methods from a number of books (including
Santamarina & Fratta, "Discrete Signals and Inverse Problems" [2005], and Menke, "Geophysical Data Analysis"  particularly chapters 3 & 4).
While the simplified forward problem exists:
(1)
d = Gm
where Okada (1985) describes the Green's function, G, that relates model parameters m, to the data, d.
If data are numerous, and/or model parameters are greater than one, each can be written as a vector of size n and m, with the G being a matrix sized n x m.
The inversion equation, where we try to solve for model parameters becomes:
(2)
Gm = d
For a leastsquares solution, the simplified inversion becomes:
(3)
m^{est} = [G^{T}G]^{1}G^{T} d
Because we expect there to be interdependence between adjacent model patches, we can describe an additional set of equations to apply to the inversion and regularize the final model.
Every model patch m_{i} can describe its model roughness as the Laplacian of that component of slip (e.g. ds, ss, ts). We will try to minimize that roughness for each slip component of the model, s, which here will be the same as m:
(4)
d_{p}= 0 = κ∇^{2}s
where κ is the smoothing kernel.
Furthermore, we combine the data weights and errors into a weighting matrix W.
Now, incorporating the errors and weighting, the forward problem becomes:
(5)
d = Gm =
$\left[\begin{array}{}\mathrm{Wiid}\\ 0\end{array}\right]$
=
$\left[\begin{array}{}\mathrm{WiiG}\\ \mathrm{\kappa \nabla 2s}\end{array}\right]$ m
where the data and greens function matrix are now everything in their respective matrices above.
The diagonal components of the weighting matrix W_{ii} = (w_{i})^{0.5}/e_{i},
which are the per datum weighting w_{i} and errors e_{i}.
With d and G now containing the smoothing and weighting terms, our model estimate can now be defined as:
(6)
m^{est} = [G^{T}G+ ε^{2}I]^{1}G^{T}d
where $\epsilon 2$I = κ∇^{2}s.
We can simplify this by defining the Generalized Inverse Matrix:
(7)
G^{g} = [G^{T}G+ ε^{2}I]^{1}G^{T}
so:
(8)
m^{est} = G^{g}d
The above mathematical formulation is modified for the bounded conditions used for the parameter estimation in lsqlin,
including the option of positivity constraints, which make the full formulation much more complicated.
However, the grouping of terms are important particularly for the development of distributed slip models as input into the solver,
and for understanding the model resolution matrix and ultimately resolution spread (an estimate of parameter importance).

We'll start with the most straightforward case of a single vertical fault plane with strikeslip motion on the interface.
Example Forward Model Run:
The simplist model includes just one fault and one datapoint ran in a forward sense to prediction deformation at that datapoint. We will not stray significantly from there.
Input:
 We will create a vertical planar fault approximately 20 km long, with 1 meter, strikeslip motion, buried between 5 and 15 km depth. The output will describe the spatial deformation every 2 km across a 60x60 km grid surrounding the fault.
 To keep the algorithm from attempting an inversion, it is important that all slip parameter ranges for a given sliptype on a fault block have equal value (i.e. ss0=ssX, ds0=dsX, and ts0=tsX in this example).
 If all range values are set to zero '0', the initial value will be used. Here, we will have 1m of leftlateral (+) strikeslip that will not vary (all range parameters are set to zero).
 Input files should all end in .in
 The input file can be created in your favorite texteditor (e.g. vim). For simplicity here, we'll just direct stdout from an echocommand into a file called local1.in. The % represents a shell prompt.
% echo "
coord local # otherwise it will expect geographic coordinates
#fault type name lon1 lat1 lon2 lat2 z1 z2 dip ss ds ts ss0 ssX ds0 dsX ts0 tsX Nd Ns
fault 2 myfault 0 10e3 0 10e3 5e3 15e3 90 1 0 0 0 0 0 0 0 0 1 1
#grid name Erot Nrot lon1(x1) lat1(y1) lon2(x2) lat2(y2) Ne Nn
grid 1kmx1km 0 0 30e3 30e3 30e3 30e3 31 31
" > local1.in
Run:
 To run GTdef, you must call up matlab (interactive and gui are not necessary).
 The processing code that oversees both forward and inverse models is GTdef and is called with just the input file name, and optionally, the number of processors you'd like it to use on your machine. By default, it will use all available processors if not defined.
# within matlab
> GTdef('local1.in',1) % will run on '1' core.
Output:
 This will create an output file called 'local1_fwd.out' that will repeat the model, and give predicted results for the defined grid locations (now as 'point' datatype).
Plotting results:
 Quick and dirty visualization can easily be achieved using the freely available gnuplot program. This allows us to look at the full 3D vector field.
# first comment out all lines that are not output data using vim or similar. (not shown)
% gnuplot
gnuplot> set ticslevel 0.0 # removes empty whitespace below plot
gnuplot> splot 'local1_fwd.out' u 4:5:(0):($7*1e5):($8*1e5):($9*1e5) w vect t "surface displacement" # just the surface results amplified to see variability
#If you're more adventurous, you can plot the fault too (this assumes the fault is vertically dipping).
gnuplot> splot "< awk '$1==\"point\"{print $0}' local1_fwd.out" u 4:5:(0):($7*1e5):($8*1e5):($9*1e5) w vect t "surface displacement",\
"< awk '$1==\"fault\"{printf \"%f %f %f\\n%f %f %f\\n%f %f %f\\n%f %f %f\\n%f %f %f\\n\", $4,$5,$8,$6,$7,$8,$6,$7,$9,$4,$5,$9,$4,$5,$8}' local1_fwd.out" w l lw 2 lt 1 t "fault"
# use the cursor to move around the data

(from 'examples/Uniformslip/Solomon')
Input:
We will need to build the model, including gathering the appropriate data and then estabishing an appropriate a priori interface. In this example, we will use the very small tsunami earthquake described in Newman et al. (2011).
 Data: Here, we pull coastal subsidence data from sheet 3 of the
excel spreadsheet.
The data, being just vertical uplift/subsidence data need to be convert into the 'point 1' format in GTdef.
Errors are estimated to be approximately 0.1 m on all data (assuming error in tidal correction and prior hightide estimates),
Because all data were measured using the same methods and likely have similar error values, we choose to give them all equal weight.
#point type name lon lat z Uv eUv weight
point 1 RendovaRendova_Harbor 157.33602 8.40359 0.0 0.15 0.10 1.0
point 1 RendovaEpata_Creek 157.30622 8.43730 0.0 0 0.10 1.0
point 1 RendovaMbaniata 157.26260 8.63325 0.0 0.70 0.10 1.0
point 1 RendovaHofofo_Pt 157.19633 8.56530 0.0 0 0.10 1.0
point 1 RendovaHabila 157.22920 8.60414 0.0 0.60 0.10 1.0
point 1 RendovaRava_Pt 157.40336 8.72264 0.0 0.60 0.10 1.0
point 1 TetepareTofa 157.53432 8.75576 0.0 0.40 0.10 1.0
point 1 TetepareJetty_near_Ecolodge 157.44286 8.72234 0.0 0.25 0.10 1.0
point 1 TetepareEcolodge_boat_ramp 157.44321 8.72120 0.0 0.30 0.10 1.0
point 1 RendovaRano 157.32886 8.62969 0.0 0.50 0.10 1.0
point 1 RendovaVankuva 157.33953 8.60934 0.0 0 0.10 1.0
point 1 RendovaKofi_Bay_village 157.33874 8.6039 0.0 0.40 0.10 1.0
point 1 RendovaMauru_Loging_Camp 157.39881 8.5137 0.0 0.30 0.10 1.0
point 1 RendovaUgele 157.39921 8.44959 0.0 0 0.10 1.0
Modeled Fault: We need to come up with a choice for the appropriate input model. This can be done using existing published fault geometries, evaluating distributions of microseismicity, focal mechanisms, examination of surface topography/bathymetry, or some combination of these methods. Once an initial model is chosen, output from a number of inversions using differing geometries can be used to chose an optimal geometry by evaluating each geometry's overall model fit. The output results of these inversions can then be broken up into distributed models (discussed below).
Here, Newman et al. (2011) defined the fault using a combination of the overall fault strike, and testing a range of dip geometries.
For our inversion, we will define a finite range of allowble thrust(+) values (ds, and dX between 0 and 100 m).
Because data here are only vertical, we will chose to set strikeslip to 0 and not allow it to vary in the model, as this component will be very poorly constrained, possibly being dependent on choice of geometry. Finally, because this is a uniform slip model, the number of slip and strike components must be set to one (Nd=Ns=1).
Note: Here, the dip =158° for a strike of 125°. This is equivelent to a dip of 22° striking 305°, but from the far endpoint.
#fault type name lon lat z1 z2 len str dip ss ds ts ss0 ssX ds0 dsX ts0 tsX Nd Ns
fault 1 slm 157.0990 8.6920 0 5200 50000 125.0 158.0 0 0.1 0 0 0 0 100 0 0 1 1
Run: After combining the above parameters into a single input file (e.g. 'solomon.in'), you can run the code just as was done in the forward model.
Output:
 Output files: If the model converges, you can analyze the output results in either the '_inv.out' or '_kp*.out'.
The former file gives the fit results of the run, while the later gives the model results as well as their predicted displacements at the sites described in the orginal input model. Its this second file that can be used as an input for future runs.
 Model fit: Below is the output from 'solomon_inv.out' created by GTdef. Since there is only one inversion ran, there will be only one data line following a header line preceded by its column number (#).
#(1)beta (2)kappa (3)data_num (4)slip_num (5)ndf (6)rss [m^2] (7)rms [m] (8)wrrs [m^2] (9)wrms [m] (10)chi2 (11)rchi2 (12)r_1d [cm/km] (13)r_2d [cm/km^2] (14)strain [cm/km]
0.00000e+00 0.00000e+00 14 1 13 4.05607e01 1.70211e01 4.05607e+01 1.70211e+00 4.05607e+01 3.12005e+00 NaN NaN NaN
 Model results and predicted displacements: Shown here is the output from 'solomon_kp0.00000.out'. Note that in addition to lines repeating the input file, there are several default input options explicitely written out here (e.g. coord, earth, surface). Also included are model misfits repeated from the above file. Lastly, the model line now includes the model output for a bestfit (showing 5.16091 m of thrust), and the model predictions for each data line (e.g. RendovaMbaniata predicts 62.62 cm of subsidence).Though data type was input as point 1, they are output as point 3, giving the predicted horizontal motions as well. Errors are not reported for models.
coord geo
#data_num 14
#slip_num 1
#ndf 13
#rss 4.05607e01 [m^2]
#rms 1.70211e01 [m]
#wrss 4.05607e+01 [m^2]
#wrms 1.70211e+00 [m]
#chi2 4.05607e+01
#rchi2 3.12005e+00
#r_1d NaN [cm/km]
#r_2d NaN [cm/km^2]
#strain NaN [cm/km]
earth homogeneous 3.00e+10 0.2500
kappa 0.00000
beta 0.00000
smooth 2d
surface free
fault 1 slm 157.09900000 8.69200000 0.0000e+00 5.2000e+03 5.0000e+04 125.00 158.00 0.00000 5.16091 0.00000 0.0000 0.0000 0.0000 100.0000 0.0000 0.0000 1 1
point 3 RendovaRendova_Harbor 157.33602000 8.40359000 0.00000e+00 0.09951 0.25301 0.00474 NaN NaN NaN 1.00
point 3 RendovaEpata_Creek 157.30622000 8.43730000 0.00000e+00 0.11693 0.31236 0.00387 NaN NaN NaN 1.00
point 3 RendovaMbaniata 157.26260000 8.63325000 0.00000e+00 0.89590 1.40681 0.62620 NaN NaN NaN 1.00
point 3 RendovaHofofo_Pt 157.19633000 8.56530000 0.00000e+00 0.26215 0.92622 0.30033 NaN NaN NaN 1.00
point 3 RendovaHabila 157.22920000 8.60414000 0.00000e+00 0.70254 1.31715 0.55581 NaN NaN NaN 1.00
point 3 RendovaRava_Pt 157.40336000 8.72264000 0.00000e+00 0.98978 1.39069 0.60162 NaN NaN NaN 1.00
point 3 TetepareTofa 157.53432000 8.75576000 0.00000e+00 0.71027 0.73619 0.19750 NaN NaN NaN 1.00
point 3 TetepareJetty_near_Ecolodge 157.44286000 8.72234000 0.00000e+00 0.87663 1.18799 0.39089 NaN NaN NaN 1.00
point 3 TetepareEcolodge_boat_ramp 157.44321000 8.72120000 0.00000e+00 0.86957 1.17774 0.38133 NaN NaN NaN 1.00
point 3 RendovaRano 157.32886000 8.62969000 0.00000e+00 0.73815 1.12756 0.30286 NaN NaN NaN 1.00
point 3 RendovaVankuva 157.33953000 8.60934000 0.00000e+00 0.59654 0.93641 0.17707 NaN NaN NaN 1.00
point 3 RendovaKofi_Bay_village 157.33874000 8.60390000 0.00000e+00 0.57048 0.90517 0.16057 NaN NaN NaN 1.00
point 3 RendovaMauru_Loging_Camp 157.39881000 8.51370000 0.00000e+00 0.24266 0.41934 0.00177 NaN NaN NaN 1.00
point 3 RendovaUgele 157.39921000 8.44959000 0.00000e+00 0.15468 0.29970 0.00941 NaN NaN NaN 1.00
 Predicting displacements elsewhere: To estimate model displacements outside of the original data points in '_kp*.out', you can copy the model output into a new file along with the points that you want prediced. Here, we are just using a grid, similar in the foward model example above, name it something like 'solomon_fwd.in', and run as before.
#fault type name lon lat z1 z2 len str dip ss ds ts ss0 ssX ds0 dsX ts0 tsX Nd Ns
fault 1 slm 157.09896 8.69204 0.00e+00 5.20e+03 5.00e+04 125.00 158.00 0.000 5.167 0.000 0.00 0.00 0.00 0.00 0.00 0.00 1 1
#grid name Erot Nrot lon1 lat1 lon2 lat2 Ne Nn
grid Solom_region 0 0 156.4 9.3 158.1 7.9 200 200
Plotting results:
 One tool that is particularly useful is the supplied algorithm 'GTdef_project', as it can take the predicted output model and transform it into a highlyreadable file for your favorite plotting program.
> % within matlab
> GTdef_project('solomon_kp0.00000.out')
This creates a new file (here 'solomon_kp0.00000_patches.out') that describes all elements (here only one) of the model, including geographic position of corners, the element center position, and slip values.
While you could use a more fullfeatured plotting program like GMT, we'll again use gnuplot for our quick and dirty plot to show the predicted displacements from '_kp*.out', and then model geometry from '_kp*_patches.out'.
# within gnuplot
# Note: We are using diffential scaling of the displacement vectors to scale with the geographic position and depth (in degrees° and km)
gnuplot> splot "< awk '$1==\"point\"{print $0}' solomon_kp0.00000.out"\
u 4:5:(0):($7*1e1):($8*1e1):($9*1e1) w vect t "surface displacement",\
"< awk '$1==\"slm\"{printf \"%f %f %f\\n%f %f %f\\n%f %f %f\\n%f %f %f\\n%f %f %f\\n\", $4,$5,$6/1000,$7,$8,$9/1000,$10,$11,$12/1000,$13,$14,$15/1000,$4,$5,$6/1000}' solomon_kp0.00000_patches.out"\
w l lw 4 lt 7 8 t "fault"
 After using your pointer to move around your model appropriately, your result should look something like this:

(from 'examples/Nonuniformslip/Solomon')
This problem is a little more tricky, as we will quickly go from an overdetermined problem (more data than model paremeters),
to an underdetermined problem.
In order to get around this, and to solve for realistic slip, we constrain the interdependence between individual slip patches.
The most popular way to do this is to define minimize the 'roughness' of the slip patch, as defined in
Harris and Segall (1987),
but more concisely in Jónsson_et al_(2002).
The method works to minimize the 2D secondorder derivative of the slip surface so that large changes are constrained within a certain value.
The way it's been implemented here is that a series of roughness solutions are determined, all yielding different misfits to the data.
It is a somewhat arbitrary decision to chose the optimal roughness parameter. Usually this is done by evaluating the kinkpoint by which
additionally smooth results significantly increase misfits.
Input
 Choose your kappa value (parameter that controls roughness) to range between 0 to 5000 by 500.
 Setup your fault to now have 70 subfaults, 7 patches alongstrike and 10 patches alongdip.
#fault type name lon lat z1 z2 len str dip ss ds ts ss0 ssX ds0 dsX ts0 tsX Nd Ns
fault 1 slm 157.0990 8.6920 0 5200 50000 305.0 22.0 0 0.1 0 0 0 0 100 0 0 7 10
Run
 Within Matlab, run GTdef('solomon_inv.in') like in the example above except this time you're running an inversion due to the way the fault parameters are set up.
Output
 This time, you have two kinds of output files, those that end in '_kp*' correspond to each kappa value you set in step "i" above, and the other is a summary file '_inv.out'. The inversion summary file will be used in the next step to determine the optimal smoothing factor (kappa value) for this dataset.
 Explore your inversion results. Using GTdef_misfit_kappa.gmt, try plotting misfit vs. roughness to see how your choices in kappa affect these results. Looking at the Lcurve, the 4th lowest kappa value is optimal for this example dataset.
# in shell (command line)
GTdef_misfit_kappa.gmt M nicoya_camp_inv_inv.out
Plotting results:
 Finally, once you've chosen your optimal kappa value, run a new forward model, and repredict deformation. The new forward model input will be a modified version of the corresponding '_kp.out' file to the optimal kappa value determined from the step above. Remember, for a forward model, it is important that all slip parameter ranges for a given sliptype on a fault block have equal value (i.e. ss0=ssX, ds0=dsX, and ts0=tsX).
 You might want to look at your results in a geographic reference system. We've created a GMT script (GTdef_slipmodel.gmt) to do this,
but you may need to modify it for your purposes.
# in shell (command line)
GTdef_slipmodel.gmt M solomon.in K750 R157.0/157.7/9/8.3 S0.1/1/0.1 D
Here is an example of a plot from the results obtained with an inversion with kappa=750.

(from 'examples/Nonuniformslip/Nicoya/Checkboard')
Develop a Synthetic Forward Model:
 Define the size boxes that you will want to try to recreate in your checkerboards.
As a general rule, we find that the model patch size that you can interpret is dependent on and about equal to its depth relative to stations overhead.
This assumes stations are immediately overhead.
This will certainly change as a function fault dip, and station coverage, thus while you may have theoretical fidelity to recover even extremely small
patches near the trench (where the fault is just below you), you may not have data until you're on land, 50 to 200 km away.
Thus, we will take an initial example of a fault patch that is 50 km laterally from a 20° dipping fault surface,
the checkerboard patch should be around 50km*sin(20°)=17km . At 100 km, it's about 34 km, and at 200 km it's about 70 km, assuming constant dip.
Note, that this is just a starting and approximate parameter.
 You will now need to define a number of subfaults that match your expected optimal resolution.
This too can be done by hand, but the program, GTdef_ckbd_subflt.m makes it much easier.
Check out the header notes within the program for details.
 For a fault named 'nicoya' with 30x40 patches and approximately 200 km x 200 km in size,
you can run
% within matlab
GTdef_ckbd_subflt('nicoya',30,40,3,4, [0 .0823 0 0 0 0 0 0 0],[0 0 0 0 0 0 0 0 0])
to define a fault that groups alternating blocks into 3 along dip, and 4 along strike patches. Here, the
the blocks are 20x20 km (smallest size for coastal datathis is too small for deeper patches).
The vectors here define the forward behavior of the patches that we normally think of as 1's and 0's.
So, patch 1 (the black, or slipping patch), p1= is only backthrust at 82.3 mm (no inversion set for these patches).
For patch 0 (the white, or null patch), p0= zeros (9 of them).
The program will create a file called GTdef_ckbd_subflt.out that will only contain the subfault parameters.
This will need to be copied into your input model file that contains your fault characterization
and 'points' in which you want to predict deformation for the forward run.
# in shell (command line)
cat GTdef_ckbd_subflt.out >> nicoya_ckbd.in
 Once it's ready, run
% within matlab
GTdef('nicoya_ckbd.in')
If everything works you will get your initial noisefree synthetics.
 Using GTdef_slipmodel.gmt, plot the noise free model with the nicoya_ckbd_fwd.out
(the output from the forward run).
# in shell (command line)
GTdef_slipmodel.gmt M nicoya_ckbd.in S100/100/10
The figure should look like the example below.
Note, that because this program is written for real data, the labeling may not be accurate. Here, 'observed' is arbitrary, and
whatever was in your original input file. Modeled data are the predicted results that will be used for subsequent inversions.
Adding Noise:
Good synthetic data includes an approximation of the noise that you'd get in the real data.
For this, We've developed a simple shell script that searches only point data from an output
file and adds Gaussian noise to the data itself.
The program then writes out the perturbed data and error estimates representing the input error,
of the size defined within the script.
This program needs updating for easier and general use, but it's pretty easy to run, as is.
 Copy the program to your working directory (its called add_noise.sh and is in the
plotting directory shared with the GTdef scripts).
 Edit the parameters at the top of the file to include the expected
East, North, and Vertical noise in your data (in meters).
Now run it it on the output file, directing into a new output.
As an example:
# in shell (command line)
add_noise.sh nicoya_ckbd_fwd.out > nicoya_ckbd_fwd_wNOISE.in
Note that the directed output is now ending in '.in'.
Invert Synthetics:
Here we will determine what the synthetic and noised data show on the fault.
 Edit new synthetic input file to either comment out or remove the subfault definitions.
You can do this manually, or run a script:
# in shell (command line)
# this scrubs file directing it into a new file that then replaces the original.
grep v subfault nicoya_ckbd_fwd_wNOISE.in > j1 ; mv j1 nicoya_ckbd_fwd_wNOISE.in
 You can add or modify any other parameters that you want into the header of the file input file.
We want to be sure that we're running a range of kappa values.
Here it's 5 between 10005000, 5 between 7k and 15k, and 4 between 20k and 50k, but the useful range
depends on the number of patches and relation to data.
Modify nicoya_ckbd_fwd_wNOISE.in to include:
surface free # slip to the trench here
kappa 2 1000 5000 5 # 1k5k for with n=5
kappa 2 7000 15000 5 # every 2k now
kappa 2 20000 50000 4 # every 10k now
resolution 1 # output the results of the resolution matrix (we can use later)
lsqlin 10 1e10 # less iterations and higher tolerance limits (faster runs with no real degradation).
project on # runs through and creates easy to read geographically projected outputs of model and slip parameters.
Change the fault line to invert (e.g. make thrust range from 0.823 to 0, the ds0 and dsX values).
#fault is set for inversion (ds0, dxX = 0, .0823)
#fault type name lon1 lat1 lon2 lat2 z1 z2 dip ss ds ts ss0 ssX ds0 dsX ts0 tsX Nd Ns
fault 2 nicoya 85.54233907 8.74026277 86.82826920 10.01077299 0 6.2e+04 17.0 0 0.0823 0 0 0 0.0823 0 0 0 30 40
 Run the inversion:
% within matlab
GTdef('nicoya_ckbd_fwd_wNOISE.in')
 Determine appropriate Κ value to use by evaluating the information in 'nicoya_ckbd_6x8_fwd_lowNOISE_inv.out', which
can be plotted with GTdef_misfit_kappa.gmt.
# in shell (command line)
GTdef_misfit_kappa.gmt M nicoya_ckbd_3x4_fwd_lowNOISE_inv.out
Results are
here.
 Now, evaluate/visualize the results of the inversion to determine which patches have better or worse representation
of the original input model that you created (your checkerboards).
A great way to do this is to create very similar plots of both the result of the forward prediction (nicoya_ckbd_fwd.out)
and the new inverse at the kappa that you feel most appropriate (e.g. nicoya_ckbd_fwd_wNOISE_k9000.out).
Now the files should be ready for the helper program GTdef_slipmodel.gmt.
# in shell (command line)
GTdef_slipmodel.gmt M nicoya_ckbd_fwd_lowNOISE.in K 9000 S100/100/10
The figure should look like the example below.
Notes:
As you can see in the above figure, the checkerboards are only replicated near the coast, and the predicted model becomes
smooth with apparent larger patches both updip and downdip, denoting loss of resolvability from a lack of stations immediately overhead (updip),
and insufficiently large patches sizes (downdip).
To improve on the downdip problem, it is possible to make a hybrid checkerboard model with larger downdip patches
(maybe 40x40 km), to test the resolution of larger slip features.
When doubling the patch sizes, features are much better resolved (in 'examples/Nonuniformslip/Nicoya/Checkboard/Nicoya6x8patches').
 input checkerboard
 tradeoff curve
 inversion results
 Checkerboard testing is relatively straightforward, but a somewhat cumbersome problem with a visual and not very quantitative
assessment of the model quality given data. Instead, it may be much better to look at the output of the resolution matrix, itself.
 The examples used in the checkerboard tests above use campaign data from the Nicoya Seismic Cycle Observatory GPS data.

Theory:
If you've ran inversions with the 'resolution' flag set to one (see option in Checkerboard Testing section), we
should have outputs with the tail '_patches_R.out' for each of the kappa values in our models. This file describes the output of
the diagonals of the model resolution matrix, R_{jj}, as the product of the Generalized Inverse Matrix and the Green's function (see Theory section).
(9)
R = G^{g}G
Equation(7) from above can be rewritten as:
(10)
G^{g} = [G^{T}W_{e}G+ ρW_{m}]^{1}G^{T}W_{e}
where W_{m} is the Laplacian of the slip, W_{e} is the error weighting matrix, and ρ is the roughness intensity applied to the models
(Jonsson et al., 2002; Kyriakopoulos and Newman, 2016).
From the diagonals of the resolution matrix, R_{jj},
you can get a characterization of how spatially a model patch's resolution is spread output across the model space.
This can be related in a single scalar per patch that is dependent on the patch size.
This parameter is called the resolution spread, r_{j}, and is
(11)
r_{j} = L_{j} / √R_{jj},
where L_{j} is the nominal patch size (Funning et al., 2005; Kyriakopoulos and Newman, 2016).
Because R doesn't depend on the actual values of the data, but only its location and weighted errors, relative to the model,
the same is true for r_{j}.
Evaluation:
(from 'examples/Nonuniformslip/Nicoya/Checkboard')
We've created a program that calculates L_{j}, and r_{j} for a selected component of slip.
The program GTdef_Rmodel.gmt uses similar inputs to GTdef_slipmodel.gmt, but requires Kappa (K) and Slip components (C) to be described.
As noted above, the program requires the respective '_patches_R.out' files that are created with the c>'resolution' flag set to 1 before an inversion.
# in shell (command line)
GTdef_Rmodel.gmt M nicoya_ckbd_fwd_lowNOISE.in K7000 S D D
This will plot the dipslip (S D) component, at Κ=7000, displaying the result (D).
The figure should look similar to the example below.
The real beauty of the Resolution Spread is in it's simplicity and quantitative nature. Instead of creating apriori assessments
of approximate resolvable model size, developing a model/running it, and doing "squintyeyed" qualitative analysis of what you can
resolve (as is done with Checkerboard testing), you get real numbers, with direct resolvable properties of patch size.
While plotting r_{j} gives a direct image of all patches average spread, you can use resolution 2 to
get information on the offdiagonal components of the model resolution matrix to plot the actual spatial smearing of individual
patches, per component (see Kyriakopoulos and Newman, 2016). You may need to develop your own plotting algorithm for this.

In the examples above, we've been working with a homogeneous elastic halfspace. In this section, we will describe another set of algorithms for creating layered models within GTdef.
The algorithms GTdef_edgrn.m and GTdef_edcmp.m, create the Green's functions and compute the expected deformation, respectively.
The codes are adopted from Wang et al.(2003), originally written in FORTRAN.
Input
(from 'examples/Layered_Models/Nicoya/All_stations')
To build an input file to run the layered model it is necessary to set the earth structure to be 'layered'
and define appropriate parameters for the desired layers.
Details are included in the example input file. The layer parameters are defined as follows;
# nl  number of lines of layer information (!= number of layers)
# each layer except the top layer (interface of halfspace) has two lines of records
# if upper & lower values are different, layers will be treated as constant gradient
# obsz  uniform depth of observation points [m]
# nr  number of radial points for green functions
# minr  minimum radial distance [m]
# maxr  maximum radial distance [m]
# nz  number of depth points for green functions
# minz  minimum depth (positive points down) [m]
# maxz  maximum depth (positive points down) [m]
# srate  sampling rate for wavenumber integration
# the ratio between the Nyquist wavenumber and the really used wavenumber
# suggested value 10128 {12}
# the larger the value is, the more accurate the results are, more time required
#earth layered nl obsz nr minr maxr nz minz maxz srate
earth layered 22 0.0 1201 0.0 300e3 326 0.0 65e3 12
Run
For this example, we use a combination of Campaign and Continuous stations from the Nicoya dataset.
The layers are defined by Vp, Vs, depth and density (rho),
in this example we use an initial 1D bestfit velocity model specific to the region as published in
DeShon et al.(2006).
In this example, each layer has internally constant properties, and are described from shallowest to deepest.
The bottom layer defines everything underneath it.
Keep in mind that each line is the layer boundary conditions, and not conditions across the entire layer.
A complete input file can be found in the examples directory (Layered_models). The first 5 lines are shown below:
#layer id depth Vp [m/s] Vs [m/s] rho [kg/m^3]
layer 1 0.0e3 5.35e3 3.01e3 2646.291
layer 2 5.0e3 6.12e3 3.44e3 2736.762
layer 3 5.0e3 6.12e3 3.44e3 2736.762
layer 4 9.0e3 6.12e3 3.44e3 2736.762
layer 5 9.0e3 6.12e3 3.44e3 2736.762
After the model is set up, the model can be run the same way as other models:
% within matlab
GTdef('nicoya_all_stations_inv.in')
Output
Since this is set up to run an inversion, the output files will be the same type ('_kp*.out')
however it too will create outputs for the algorithms' Green's functions by default.
The GTdef_edgrn.m algorithm will generate these and store them in a subdirectory called ./edgrnsfcts.
The directory will contain three large files describing each strikeslip (ss), dipslip (ds), and the compensated linear vector dipole (cl).
Once the files are created, you can change the fault and station geometry minding not to go beyond the depth and distance limits established
in the earth layered line.
Plotting results
Using the gmt plotting programs introduced in the sections above,
we can plot the misfit curve and the geographic representation of the slip and deformation estimates.
# in shell (command line)
GTdef_slipmodel.gmt M nicoya_allstations_layered_inv.in K15000 D
For the Nicoya Dataset and the layered model, a representative result is shown below.
For comparison, we include the Okada solution
for the same fault geometry and data.

(from 'examples/Joint_Inversion/3_Tsunami')
(to be complete)

 Santamarina, J. C., & Fratta, D. (2005). Discrete signals and inverse problems. An Introduction for Engineers and Scientists. UK: Wiley & Sons.
 Menke, W. (2012). Describing inverse problems.
 Wang, R., Martı́n, F. L., & Roth, F. (2003). Computation of deformation induced by earthquakes in a multilayered elastic crust—FORTRAN programs EDGRN/EDCMP. Computers & Geosciences, 29(2), 195207.
 Funning, G. J., Parsons, B., Wright, T. J., Jackson, J. A., & Fielding, E. J. (2005). Surface displacements and source parameters of the 2003 Bam (Iran) earthquake from Envisat advanced synthetic aperture radar imagery. Journal of Geophysical Research: Solid Earth, 110(B9).
 Kyriakopoulos, C., & Newman, A. V. (2016). Structural asperity focusing locking and earthquake slip along the Nicoya megathrust, Costa Rica. Journal of Geophysical Research: Solid Earth, 121(7), 54615476.
 Okada, Y. (1985). Surface deformation due to shear and tensile faults in a halfspace, Bull. Seism. Soc. Am. 75, 11351154.
 Chen, T., A. V. Newman, L, Feng, H. M. Fritz (2009), Slip Distribution from the 1 April 2007 Solomon Islands Earthquake: A Unique Image of NearTrench Rupture, Geophys. Res. Lett., 36, L16307, doi:10.1029/2009GL039496.
 Feng, L., A. V. Newman, J. M. Protti, V. González, Y. Jiang, T. H. Dixon (2012), Active Deformation near the Nicoya Peninsula, Northwestern Costa Rica, Between 1996 and 2010: Interseismic Megathrust Coupling, J. Geophys. Res., 117, B06407, doi:10.1029/2012JB009230.
 DeShon, Heather R., et al. "Seismogenic zone structure beneath the Nicoya Peninsula, Costa Rica, from threedimensional local earthquake Pand Swave tomography." Geophysical Journal International 164.1 (2006): 109124.
