DisModel is a set of matlab programs and a few fortran components to invert various types of geodetic data for the geometry and slip on faults in an elastic half-space. Of course, one can also run forward models and evaluate their misfit in a trial-and-error mode. The code originated at Stanford, when Paul Segall and I put together new routines to model the Loma Prieta postseismic deformation, but there are now a large number of different versions and add-ons based on con- tributions by Sue Owen, Peter Cervelli, Evelyn Price, David Schmidt, etc. etc. Evelyn Price in particular added some nice things allowing for slip-distributed models using positivity constraints. I put DisModel files on anonymous ftp at quake.geo.berkeley.edu cd outgoing/burgmann/DISMODEL here you also find example data files and GMT plotting routines (.map file is for GMT plotting of the model geometry and vectors). To keep the core programs as default names, set: setenv MATLABPATH /data/foxnsox/burgmann/matlab/matools The matools directory in DisModel.tar has the core modeling programs, mostly developed with Paul Segall back when I was at Stanford. The basic set of matlab code in /matools includes some compiled fortran code including the actual dislocation code after Okada (1985). DUZCE_BVLS and DUZCE_INSAR contain all the programs and data I used to model the 1999 Duzce earthquake GPS and InSAR data with geometry and distributed slip inversion. You want to run all the modeling in the other directory area where you may find different versions of some of the same programs as in matools, but updated for some reason. So, if you run matlab there, these will be used rather than those in the default matools area. Essentially, once you've set up your initial fault geometry file and data files (GPS displacements) you are ready to go. A few important files: FaultGeom - has the a priori fault geometry in it, as many faults as you like using: width depth dip st_lat st_lon en_lat en_lon s-slip d-slip open km km deg deg deg deg deg x x km 4000 41 0 22.000 120.807 24.873 122.000 0 0 0 16.55 16.55 89.99 40.70 29.40 40.74 30.52 -2.76 0 0 (one line per fault) FaultBounds - prescribes bounds on geometry parameters for the non-linear inversions for geometry. len, wid, depth, dip, strike, latitude, longitude 15.0000 10.0000 15.0000 35.0000 66.0000 22.8000 70.0000 30.0000 35.0000 45.0000 75.0000 136.1000 23.5000 70.6000 VolGeom contains a priori parameters for a Mogi volume source that can be used for modeling volcanic intrusions and other volume-change processes. temp.gmtlin is a line file for plotting the fault geometry in GMT using psxy: psxy temp.gmtlin -J${projection} -R${range} -P -M -K -O -W3/000/000/000 -V >> $ {psfile} In Geometry inversion, start.gmtlin is created which has the a priori geometry and temp.gmtlin will be the final optimized geometry. volums.gmtlin is for plotting location of volume source in GMT Primary program: DisModel.m - DisModel_smooth.m - Type "DisModel" in the matlab command line and the GUI to run the whole code comes up giving you the means to load in data, set up the geometry, and run inversions. For some stuff, it's easier to run the code on command line ForModel.m - Does Forward model, using the prescribed fault geometry and slip in ForModel dhat = G*slip; r = W*(d-dhat); what is G, d W.? d is the vector of data (e.g., GPS measured displacement components W is the weight matrix, directly derived from the data errors (cov, covariance matrix). G is the matrix of Greens functions created by MakeGgps.m that relates unit slip on each fault patch to the displacements at each site computed with dislocation algorithm. slip contains the slip components (strike-slip, dip slip, opening) that are either read from FaultGeom (ForModel) or estimated via least squares (SlipEst) Thus, dhat = G*slip is the vector of predicted surface displacements. r represents the sum of the residuals (data d - model prediction dhat) SlipEst.m - Does least squares inversion for slip on all dislocations (only for those slip values that have non-zero a priori estimates). For the Izmit-EQ study I wrote a new version SlipEst_smooth.m, which allows us to smooth the distributed slip models. Evelyn Price improved on that and created the version with positivity constraints I used in DUZCE_INSAR. GeomEst.m - Non-linear inversion for fault geometry with bounds on parameters and constraints if needed. (bounds get defined in FaultBounds) The matlab optimization program constr.m is not supported anymore, so it would be advisable to use fmincon.m instead, which is implemented in the second series of programs, but I have not tested it myself. LoadData.m organizes the data (Calls Get???data scripts and builds the data and covariance matrices). GetGMTdata.m - Reads 2D GPS velocities only, from simple GMT format. (2D velocities in GMT format, but no site names in last column). GetGPSdata.m - Reads 3D velocities/displacements and full covariance. GetSARdata.m - Reads the subsampled InSAR data. MakeGSAR.m/MakeGgps.m sets up the G matrix relating unit slip/opening on the dislocation and surface displacements or range change. Of course, we can also just feed it Green's functions provided by other code, such as the slip- to-site displacement relations computed by Liz Hearn from her FEM models. The model displacements can be visually compared with the observed displacements by executing a gmt script (e.g., duzce.map) which will plot the obs.gmtvec and mod.gmtvec displacement files that we create in the model programs. (can also plot vertical (mod_vert.gmtvec etc.) To do a test run for a simple inversion using a single dislocation plane, you do the following (the file names might not be quite up to date) To run the GUI, So, to run things, (1) go into the working directory. Edit the fault geometry file, data file name in DisModel.m, and create a data file from your GPS velocities. Usually, the last station is used as a reference site and its velocity subtracted from all the others. (2) start matlab (3) run DisModel (4) click button GPS data in lower left to select data files and type the correct file name and Hit Load Data button (5) Hit Load Data button (6) Hit Setup Geometry button (the fault parameters that appear can be edited right here and saved if you changed them, useful as you can make sure your fault elements are correctly located). (7) Do a forward model or slip inversion (Hit Estimate Slip) (8) Update the GMT file for your local area etc. and run to plot the model and GPS data, as well as the surface trace of the model fault (temp.gmtlin). The plot includes both horizontal (arrows) and vertical (vertical bars) displacements if you have any included. (modvert.gmtvec and obsvert.gmtvec). Various versions of some of these programs, as well as add-ons allowing for new data types, inversion approaches etc. are in /data/foxnsox/burgmann/matlab