Example 4: Static finite-fault estimation

In this example we will determine a variable slip distribution for the L’aquila 2009 earthquake by using static InSAR data. The data is the exact same from Example 3, where the overall geometry of the fault plane was estimated. It is a requirement to have Example 3 completed in order to follow the instructions and commands given in this example.

Please make sure that you are one level above the Laquila project folder (created earlier).:

cd $beat_models_path

Init

In this example you will first make use of the mode argument. The default is: geometry, which is why it was not necessary to specify it in the earlier examples. Now we will always have to set mode to ffi, which is an abbreviation for finite-fault-optimization. The following command will create a configuration file for the ffi mode called config_ffi.yaml right next to the config_geometry.yaml:

beat init Laquila --mode='ffi' --datatypes=geodetic

It will load the config_geometry.yaml and port arguments that have been specified before to ensure consistency and will only use geodetic data.

The main differences in the two configuration files are in the geodetic_config.gf_config and the problem_config. You may want to have a first glance at the new config. The general structure is the same and the argument names are chosen to give the user an initial idea what these are for. A short explanation for each argument is again given in the module API here, where you can use the search function of your browser to find the argument of interest.

Calculate Greens Functions

For the distributed slip optimization a reference fault has to be defined that determines the overall geometry. Once this has been done the problem becomes linear as the only unknown parameters are the slips in rake perpendicular and rake parallel direction. The fault geometry needs to be defined in the geodetic.gf_config.reference_sources.:

gf_config: !beat.GeodeticLinearGFConfig
  store_superdir: /home/vasyurhm/GF/Laquila
  reference_model_idx: 0
  n_variations: [0, 1]
  earth_model_name: ak135-f-average.m
  nworkers: 4
  reference_sources:
  - !beat.sources.RectangularSource
    lat: 42.29
    lon: 13.35
    north_shift: 5542.073672733207
    east_shift: 10698.176839272524
    elevation: 0.0
    depth: 2926.702988951863
    time: 2009-04-06 01:32:49.19
    stf: !pf.HalfSinusoidSTF
      duration: 0.0
      anchor: 0.0
    stf_mode: post
    strike: 144.48588106798306
    dip: 54.804317242125464
    rake: -114.58259929068664
    length: 12219.56636799338
    width: 9347.802691161543
    velocity: 3500.0
    slip: 0.5756726498300268
  patch_widths: [2.0]
  patch_lengths: [2.0]
  extension_widths: [0.6]
  extension_lengths: [0.4]
  sample_rate: 1.1574074074074073e-05

The values shown above are parts of the MAP solution from the optimization from Example 3. The results can been imported through the import command specifiying the –results option. We want to import the results from the Laquila project_directory from an optimization in geometry mode and we want to update the geodetic part of the config_ffi.yaml:

beat import Laquila --results=Laquila --mode='geometry' --datatypes=geodetic

Of course, these values could be edited manually to whatever the user deems reasonable.

Warning

The reference point(s) on the reference_fault(s) are the top, central point of the fault(s)! Ergo the depth parameter(s) relate(s) to the top edge(s) of the fault(s).

Now we need to specify the dimensions of the patches we want to discretize the reference fault(s) into. The patch_widths and patch_lengths arguments have to be lists and the respective entry needs to result in square patches. The parameters extension_withs and extension_lengths specify by how much the reference fault(s)should be extended in each direction. Example: 0.1 means that the fault is extended by 10% of its width/length value in each direction and 0. means no extension.

Once we decided for the discretization and the reference fault values we can create the discretized fault through.:

beat build_gfs Laquila --datatypes=geodetic --mode=ffi

This will create a directory ($linear_gfs): Laquila/ffi/linear_gfs where the fault geometry is saved as a pickle file ‘fault_geometry.pkl’. We can inspect the geometry of the resulting extended discretized fault wrt. the reference fault with:

beat check Laquila --what=geometry --mode=ffi

This will open an interactive 3d plot of the fault geometry, which looks along the lines of

../_images/Laquila_FaultGeometry1.png

The grey rectangle shows the geometry of the fault specified under reference_sources and the red rectangle(s) show the extended fault with the respective discretization of the sub-patches. The grey and red dots mark the centres of the reference_fault(s) and the extended faults, respectively. The numbers are the indexes of the repsective sub-patch in the Green’s Function matrix we are going to calculate next.

Note

If the upper edge of the fault would intersect the surface (no topography assumed) due to the extension it is truncated at the intersection and not extended further. Which is why the extent of the red fault is assymetric around the grey reference fault in dip-direction.

To repeat the fault discretization after changing some parameters please add the –force option and the previously fault geometry will be overwritten.:

beat build_gfs Laquila --datatypes=geodetic --mode=ffi --force

The next command starts the calculation of the linear Green’s Function matrixes (also called library) using nworkers CPUs in parallel with unit slip in each slip-direction.:

beat build_gfs Laquila --datatypes=geodetic --mode=ffi --execute

Note

The slip components are not dip-slip and strike-slip, but rake-parallel (uparr in config_ffi.yaml priors) and rake-perpendicular (uperp in config_ffi.yaml priors) wrt to reference_fault(s) rake angle(s). This is following the convention of [Minson2013].

This will create two files for each GF library in the $linear_gfs directory:
  • geodetic_uparr_static_0.traces.npy a numpy array containing the linear GFs
  • geodetic_uparr_static_0.yaml a yaml file with the meta information

Now we are ready to prepare the optimization setup.

Optimization setup

Under the problem_config we find the parameters that we need to adjust:

problem_config: !beat.ProblemConfig
  mode: ffi
  mode_config: !beat.FFIConfig
    regularization: none
    npatches: 121
    initialization: random
  source_type: RectangularSource
  stf_type: HalfSinusoid
  decimation_factors:
    geodetic: 1
    seismic: 1
  n_sources: 1
  datatypes: [geodetic, seismic]
  hyperparameters:
    h_SAR: !beat.heart.Parameter
      name: h_SAR
      form: Uniform
      lower: [-20.0]
      upper: [20.0]
      testvalue: [0.0]
  priors:
    uparr: !beat.heart.Parameter
      name: uparr
      form: Uniform
      lower: [-0.05]
      upper: [6.0]
      testvalue: [1.15]
    uperp: !beat.heart.Parameter
      name: uperp
      form: Uniform
      lower: [-0.3]
      upper: [4.0]
      testvalue: [0.5]
  hierarchicals:
    Laquila_ascxn_offset: !beat.heart.Parameter
      name: Laquila_ascxn_offset
      form: Uniform
      lower: [-0.004496268249748271]
      upper: [-0.004496268249748271]
      testvalue: [-0.004496268249748271]
    Laquila_ascxn_ramp: !beat.heart.Parameter
      name: Laquila_ascxn_ramp
      form: Uniform
      lower: [-0.00043773457168120667, -0.00023808150002277328]
      upper: [-0.00043773457168120667, -0.00023808150002277328]
      testvalue: [-0.00043773457168120667, -0.00023808150002277328]
    Laquila_dscxn_offset: !beat.heart.Parameter
      name: Laquila_dscxn_offset
      form: Uniform
      lower: [-0.003754963750062188]
      upper: [-0.003754963750062188]
      testvalue: [-0.003754963750062188]
    Laquila_dscxn_ramp: !beat.heart.Parameter
      name: Laquila_dscxn_ramp
      form: Uniform
      lower: [4.978325480108451e-05, -0.00025072248953317104]
      upper: [4.978325480108451e-05, -0.00025072248953317104]
      testvalue: [4.978325480108451e-05, -0.00025072248953317104]

Note

The npatches parameter should not be manually adjusted. It is automatically set by running the fault discretizeation step during GF calculation(above).

Hierarchicals

Please notice the hierarchicals parameters! These are the MAP parameters for the orbital ramps for each radar scene that have been optimized in Example 2. These parameters are imported if the fit_plane parameter in the geodetic_config was set to True. The default is to fix these ramp parameters during the static distributed slip optimization, because leaving them open often results in tradeoffs with patches at greater depth and thus artificial slip is optimized at greater depth. Nevertheless, the user may want to try out to free the upper and lower bounds again to include the parameters into the optimization.

Priors

The upper and lower bounds of the two prior variables can be adjusted to reduce the solution space (slip parameters [m]). For the L’aquila earthquake it is highly unlikely to have 6 meters of slip, which is simply the default parameter. A maximum slip of 2 meters in slip parallel direction may be more reasonable. In order to be able to sample the zero value at the lower bound it is necessary to allow for some backslip- ergo negative uparr; here 0.1 might be a reasonable choice.

To also allow for variable rake angles accross the fault we may want to allow some rake perpendicular slip. Here the lower and upper bounds should be set to -1. and 1., respectively.

Note

In order to fix a variable at a certain value, the lower and upper bounds as well as the testvlue need to be set to the same value.

Regularization

The regularization argument should be set to laplacian to introduce a smoothing constraint that penalizes high slip gradients between neighboring patches. Once this is enabled we need to update the configuration file to initialize the slip-smoothing weight as a random variable in the optimization [Fukuda2008]. Adding the –diff option will display the changes to the config to screen instead of applying them to the file.:

beat update Laquila --mode=ffi --diff --parameters=hypers

Once happy with the displayed changes the changes will be applied to the file with:

beat update Laquila --mode=ffi --parameters=hypers

Note

The None regularization would be used if covariance matrices that describe the theory errors for the velocity model and/or the fault geometry have been estimated [Duputel2014] , [Ragon2018]. How to do that in BEAT will be part of another tutorial in the future.

Sample the solution space

Please refer to the ‘Sample the solution space section’ of example 3 example for a more detailed description of the sampling and associated parameters.

Firstly, we only optimize for the noise scaling or hyperparameters (HPs) including the laplacian smoothing weight:

beat sample Laquila --hypers --mode=ffi

Checking the $project_directory/config_ffi.yaml, the hyperparameter bounds show something like:

hyperparameters:
h_SAR: !beat.heart.Parameter
  name: h_SAR
  form: Uniform
  lower: [-1.0]
  upper: [5.0]
  testvalue: [2.0]
h_laplacian: !beat.heart.Parameter
  name: h_laplacian
  form: Uniform
  lower: [-5.0]
  upper: [5.0]
  testvalue: [0.5]

Markov Chain initialization

The initialization argument determines at which point in the solution space to initialize the Markov Chains. The default value random simply draws a random point in the solution space from the prior distributions for each Markov Chain to be sampled. However, as we are using a laplacian smoothing constraint we can use the non-negative least-squares solution as a starting value for a randomly drawn smoothing weight (from the initial guess on the h_laplacian parameter range)[Fukuda2008]_.

The ‘n_jobs’ number should be set to as many CPUs as the user can spare under the sampler_config. The number of sampled MarkovChains and the number of steps for each chain of the SMC sampler should be set to high values as we are optimizing now for ca 250 random variables (if the values from the tutorial haven’t been altered by the user); for example to 5000 and 400, respectively.

Warning

With these sampler parameters a huge amount of samples are going to be stored to disk! With the values from the tutorial approximately 140GB of samples are created in the course of the sampling. Please see example 0 for an instruction on how to keep only the important samples to reduce the disk usage.

Finally, we are set to run the full optimization for the static slip-distribution with:

beat sample Laquila --mode=ffi

Summarize and plotting

After the sampling successfully finished, the final stage results have to be summarized with:

beat summarize Laquila --stage_number=-1 --mode=ffi

After that several figures illustrating the results can be created.

For the slip-distribution please run:

beat plot Laquila slip_distribution --mode=ffi
../_images/Laquila_static_slip_dist_-1_max.png

To get histograms for the laplacian smoothing, the noise scalings and the posterior likelihood please run:

beat plot LaquilaJointPonlyUPDATE_wide stage_posteriors --stage_number=-1 --mode=ffi --varnames=h_laplacian,h_SAR,like
../_images/stage_-1_max1.png

For a comparison between data, synthetic displacements and residuals for the two InSAR tracks in a local coordinate system please run:

beat plot Laquila scene_fits --mode=ffi
../_images/scenes_-1_max_local_01.png

The plot should show something like this. Here the residuals are displayed with an individual color scale according to their minimum and maximum values.

For a plot using the global geographic coordinate system where the residuals have the same color bar as data and synthetics please run:

beat plot Laquila scene_fits --mode=ffi --plot_projection=latlon
../_images/scenes_-1_max_latlon_01.png

References

[Duputel2014]Duputel, Z., Agram, P. S., Simons, M., Minson, S. E., and Beck, J. L. (2014). Accounting for prediction uncertainty when inferring subsurface fault slip. Geophysical Journal International, 197(1):464–482
[Fukuda2008]Fukuda, J. and Johnson, K. M. (2008). A fully Bayesian inversion for spatial distribution of fault slip with objective smoothing. Bulletin of the Seismological Society of America, 98(3):1128–1146
[Ragon2018]Ragon, T., Sladen, A., Simons, M. Accounting for uncertain fault geometry in earthquake source inversions – I: theory and simplified application, Geophysical Journal International, 214(2):1174–1190