Dislocation Modeling using GTdef for Modern Geodetic Methods

by Andrew Newman

GTdef is an implementation of the Okada (1985) dislocation model written by Ting Chen and Lujia Feng, with codes modified from disl written by Peter Cervelli. Okada (1985) can describe the surface deformation associated with rectangular dislocations of arbitrary orientation within a purely-elastic half-space. Dislocations can be described by either dip-slip, strike-slip, 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 here 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. Here, we'll step through first some forward models, and eventually inversions.



    Code:
  • Compressed Tar file of the matlab codes used here. If you use this code for your research, please reference, the above Chen et al (2009) paper. I will try to add updated example information to the codes over time, but for the time being, the below information should be sufficiently helpful.

    Starting Notes:

    • Work in Progress: Again, we'll be modifying this as we move forward in order to develop a usable cookbook for future use.
    • Matlab: Algorithms are written in Matlab and require the 'Parameter estimation' toolbox for inversions.
    • Matlab Path: The programs currently work on our internal Linus network. Be sure to add the program path to your start-up file for Matlab to ensure you can use the program in the future.
       % echo "path('/home/lfeng/DEF/GTdef',path);" >> ~/matlab/startup.m 
    • Current Version: Right now we are working with GTdef version 1.5. Lujia is currently working on an implementation that will allow for layered elastic geometries.

    Table of Contents:

    1. Input File: This is the only documentation other than within the code, and this webpage (which is in development). As you will be able to see there are a number of ways to both describe the input 'fault' models, as well as output data.
    2. Simple Planar Forward Model: Create A vertical planar fault approximately 20 km long, with 1 meter, strike-slip motion, and buried between 5 and 15 km depth. For the output, describe the spatial deformation every 2 km across a 60x60 km grid surrounding the deformation.
       % echo "
         coord local
         fault 2 myfault 0 -10e3 0 10e3 5e3 15e3 90 1 0 0  0 0  0 0  0 0
         grid 1kmx1km 0 0 -30e3 -30e3 30e3 30e3 31 31 " > local1.in
       

      Run GTdef
       # within matlab
       > GTdef('local1.in',1) % the '1' is for one core  
      

      This will create an output file called 'local1_fwd.out'
      Now visualize the output using gnu plot to show the full vector field.
       # first comment out all lines that are not output data using vim or similar. (not shown)
       % gnuplot
         gnuplot>    splot 'local1_fwd.out' u 4:5:(0):($7*1e5):($8*1e5):($9*1e5) w vect 
        % use the cursor to move around the data
      

    3. Inversions: GTdef performs a linear least-squares inversion of Green's functions to define the best-fitting set of models to describe a dataset with errors. The program uses the parameter optimization toolbox program called lsqlin. Thus, this program will work on GT machines right now, but likely will not run on your personal machine unless you've specifically installed the parameter optimization toolbox.
      You can read more about inversion methods from a number of books. One option is: Santamarina & Fratta (2005) in our paper repository.
      Simply enough, for an over-determined problem:
      D = g(M)+e
      where Okada (1985) describes the Green's function, g, that linearly relates the model vector M, to the data vector, D, with errors, e.


      Uniform slip model:
      1. From Newman et al. (2010) pull the coastal subsidence data from sheet 3 of the excel spreadsheet.
      2. Convert this data to the GTdef input format for vertical only data. Errors are approximately 0.1 m on all data, and all data were measured using the same methods, and likely have similar error values, thus they should all be given equal weight.
        # point         
        point 1  Rendova-Rendova_Harbor         157.33602       -8.40359   0.0  -0.15   1.0
        point 1  Rendova-Epata_Creek            157.30622       -8.43730   0.0  0       1.0
        point 1  Rendova-Mbaniata               157.26260       -8.63325   0.0  -0.70   1.0
        point 1  Rendova-Hofofo_Pt              157.19633       -8.56530   0.0  0       1.0
        point 1  Rendova-Habila                 157.22920       -8.60414   0.0  -0.60   1.0
        point 1  Rendova-Rava_Pt                157.40336       -8.72264   0.0  -0.60   1.0
        point 1  Tetepare-Tofa                  157.53432       -8.75576   0.0  -0.40   1.0
        point 1  Tetepare-Jetty_near_Ecolodge   157.44286       -8.72234   0.0  -0.25   1.0
        point 1  Tetepare-Ecolodge_boat_ramp    157.44321       -8.72120   0.0  -0.30   1.0
        point 1  Rendova-Rano                   157.32886       -8.62969   0.0  -0.50   1.0
        point 1  Rendova-Vankuva                157.33953       -8.60934   0.0  0       1.0
        point 1  Rendova-Kofi_Bay_village       157.33874       -8.6039    0.0  -0.40   1.0
        point 1  Rendova-Mauru_Loging_Camp      157.39881       -8.5137    0.0  -0.30   1.0
        point 1  Rendova-Ugele                  157.39921       -8.44959   0.0  0       1.0
        		  
      3. Build a model for a single fault with no subfaults. Use your favorite geographic viewer (geomapapp, google earth, google maps) to define the boundaries of a megathrust fault. Use the CMT focal mechanism to define the dip (If you haven't used this before, it is a good opportunity to become familiar with my code searchCMT.
      4. Define the fault to allow for a reasonable range of slip values. Note that since we only have limited vertical data here, it would be unwise to consider any transverse (strike-slip) motion on the fault. As well, given nature abhors a void, it would be unwise to choose opening on this model. Once you've set up your data and your fault with a slip range, you're ready to model.
        # fault
        #fault type name  lon       lat     z1  z2   len     str    dip     ss   ds   ts     ss0  ssX    ds0  dsX   ts0  tsX 
         fault  1   slm   157.0990  -8.6920  0  5200  50000  125.0  158.0   0    0.1   0     0    0      0    100    0   0
        # point         
        point 1  Rendova-Rendova_Harbor         157.33602       -8.40359   0.0  -0.15   1.0
        point 1  Rendova-Epata_Creek            157.30622       -8.43730   0.0  0       1.0
        point 1  Rendova-Mbaniata               157.26260       -8.63325   0.0  -0.70   1.0
        point 1  Rendova-Hofofo_Pt              157.19633       -8.56530   0.0  0       1.0
        point 1  Rendova-Habila                 157.22920       -8.60414   0.0  -0.60   1.0
        point 1  Rendova-Rava_Pt                157.40336       -8.72264   0.0  -0.60   1.0
        point 1  Tetepare-Tofa                  157.53432       -8.75576   0.0  -0.40   1.0
        point 1  Tetepare-Jetty_near_Ecolodge   157.44286       -8.72234   0.0  -0.25   1.0
        point 1  Tetepare-Ecolodge_boat_ramp    157.44321       -8.72120   0.0  -0.30   1.0
        point 1  Rendova-Rano                   157.32886       -8.62969   0.0  -0.50   1.0
        point 1  Rendova-Vankuva                157.33953       -8.60934   0.0  0       1.0
        point 1  Rendova-Kofi_Bay_village       157.33874       -8.6039    0.0  -0.40   1.0
        point 1  Rendova-Mauru_Loging_Camp      157.39881       -8.5137    0.0  -0.30   1.0
        point 1  Rendova-Ugele                  157.39921       -8.44959   0.0  0       1.0
        		  
      5. If the model converges, you can now use the output model solution as the input for a new forward model to repredict deformation for the entire region. Similar to what is done in the grid example above.
        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  100.00  0.00  0.00
        #grid name           Erot  Nrot  lon1    lat1    lon2     lat2     Ne  Nn
        #grid  Solom_region   35    35    156.0   -9.5    159.0    -8.0     200   200 
        grid  Solom_region   0    0    156.4   -9.3    158.1    -7.9     200   200 
        		


      Non-Uniform slip model:
      1. 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 second-order derviative 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 yeilding different misfits to the data. It is a somewhat arbitrary decision to chose the optimal roughness parameter. Usually this is done by evaluating the kink-point by which additionally smooth results significantly increase misfits.
      2. Choose your kappa value (parameter that controls roughness) to range between 0 to 5000 by 500.
        kappa 2 0 5000 11 
      3. Setup your fault to now have 12 subfaults, 4 patches along-strike and 3 patches along-dip.
      4. Invert, and explore your results. Try plotting misfit vs. roughness to see how your choices in kappa affect these results.
      5. Once you've chosen your optimal kappa value, run a new forward model, and repredict deformation.
      6. Finally, you might want to look at your results in a geographic reference system. I've created a GMT script to do this, but you must modify it for your purposes. As well, I haven't gotten the slip models to plot on top of the geographic ones, so you may want to do this by hand using your favorite vector graphics program like Adobe Illustrator. If you do get your version to correctly plot on the map without any one-time tricks, let me know and I will update the online version.

        The code: slipmodel.gmt
      7. Another code that may be of use is GTdef_interface.gmt which can be used to plot an individual slip model along with the mistfit/roughness trade-off. You again may need to modify this for your needs.
        % GTdef_interface.gmt *kp*.out
        % # view your results 
        % gs *.ps   # page through your numerous ps files to observe the variance with increasing smoothness
        		  


      Checkerboard Testing:

      Okay, you've gotten some solutions by now, however how are you certain that your results are of value? This is a big question that is quite difficult to rightfully answer. One quote that often comes to mind here is "All models are wrong, but some are useful" attributed to George Box, a statistician.
      In order to valididate a model, it would be useful if you had independent information, like mapped fault offsets, or tsunami results. However, in the absense of this, one must work with the data we used for the inversion in the first place. This is normally done by creating a checkerboard model and determining whether the predicted 'data' are sufficient to recreate the model.
      1. Create an input model with subfaults that are approximately as large as characteristic changes you observe in your solutions. Each subfault alternates between 0 and 1 unit of slip.
      2. Run a forward model, and predict the deformation at each of the sites you used for your original inversion.
      3. If you have error information add these values so they fractionallly scale with the difference between your original data and their error. (e.g 10 cm error for 50 cm original displacment) should now scale to 1 cm error for predicted 5 cm displacemnt).
      4. Use these synthetic data and errors to run a new inversion and compare the output to the original checkerboard. If the input model matches the inversion, then you have great resolution. Where you start to signficantly deviate from your checkerboard, you lose resolution. One could be more quantitative than this, but it is usually best to trust your judgement here--it's kept you alive this long.


Course Home | anewmangatech.edu | Updated: Thu Mar 13 16:57:49 EDT 2014