Setting up a seafloor motion forward problem

Assuming you already have a rupture model in the form of a .rupt file you prepared yourself or a .inv file from an inversion your goal is to generate displacements on a grid of points on the seafloor. These can be instantaneous or not and in the end you want to create a .dtopo dile to be used with GeoClaw

Generate grid of points

The number of grid points is restricted by the header file /src/fk/model.h by the line parameter(nlay=20, ndis=4000, nt=4096). In this case ndis=4000 means 4000 points is the maximum allowed. You can change this but then will need to re-compile fk.

Anyway, then make a .sta file with the grid points by doing:

>>> from mudpy.clawtools import make_grid
>>> make_grid(-76,-71.5,-40,-31,0.1,0.1,'TG','/Users/dmelgar/Slip_inv/maule_gps/data/station_info/seafloorGF.sta')

This will make a .sta file with stations named TG0000, TG0001, etc. between (-76,-40) and (-71.5,-31) at 0.1 degree intervals in latitude and longitude and save it to the specified path

Generate GFs and synthetics

Now do some text editing or scripting to change the .sta file from the previous step into .gflist file (sorry no automation for this yet). Load up template.fwd.py from the run/ folder or a previous .fwd.py parameter file. Set the flags

make_green=1
make_synthetics=1 
solve=1 

And make sure rupt_file and GF_list are pointing to the right files. Run the parameter file and you will obtain statics or displacement waveforms at all the grid points. Depending on the complexity of the rupture file and the number of grid points and type of output requested this can take a good long while O(10ish hrs).

When this is done you should see output in your project folder under \output\forward_models\

Generate bathymetry derivatives for topography effect

First you need to prepare your topo/bathy file. For example I have a netCDF .grd file named maule_30.grd that covers a region larger than what I need so first I cut it with GMT

$ gmt grdcut maule_30.grd -R285.75/288.45/-38.625/-33 -Gmaule_geoclaw.grd

Now make the x derivative with

$ gmt grdgradient maule_geoclaw.grd -Gmaule_geoclaw_dx.grd -A90 -N

The -A90 flag indicated the derivative is taken in the due east direction and the -N indicates it is normalized so that the amplitude is in the range +/- 1. Do the same for the y derivative (now in the north direction) by

$ gmt grdgradient maule_geoclaw.grd -Gmaule_geoclaw_dy.grd -A0 -N

You can quickly verify by plotting one of the output files using mudpy.view.plot_grd as

>>> from mudpy import view
>>> view.plot_grd(u'/Users/dmelgar/DEMs/maule/maule_geoclaw_dx.grd',[-0.5,0.5],cmap=cm.seismic)

Make a coast_path file

The topography effect should only be applied to points below the sea surface so it's important to create a path file that is very close to the coast in order to filter out which points will be affected by the topo effect and which won't. Make a path in Google Earth, save it as a .kml not a .kmz and then convert it to text with GMT by doing

$ gmt kml2gmt coast_path.kml > coast_path.txt

Make sure you edit out the > character in the header that GMT creates

Make .dtopo file

Now we need to process the grid output into a .dtopo file that GeoClaw can parse. For this purpose we will use mudpy.forward.move_seafloor, this function needs the topo/bathy x and y derivative grids as input in order to apply the horizontal advection effect. Make sure you point topo_dx_file and topo_dy_file to the correct .grd absolute paths. Also point coast_file to the coast path file made with Google Earth. tsun_dt is the sampling rate at which you created the seafloor synthetics and maxt the maximum time you want output at. Once all these variables are defined simply run

from mudpy.forward import move_seafloor
move_seafloor(home,project_name,run_name,topo_dx_file,topo_dy_file,tgf_file,
        outname,time_epi,tsun_dt,maxt,coast_file,topo_effect=True,variance=None,static=False)

and that ought to do it. If the area is large and the time series long this can make a rather big file (several GB) so be weary. When the calculation is done you should see a .dtopo file in your project folder under /output/forward_models/.