Preliminary eigenmode approach working with zero slip boundary conditions

[brendanjmeade]

@jploveless please take a look at https://github.com/jploveless/relax_disloc/blob/main/relax_disloc_eigen_no_zero.ipynb.

Things that this does:

• Classical smoothed solution (without zero slip BCs)
• Eigenmode solution (with zero slip BCs)

Things that I plan to do:

• Clean this up to make a nice demonstration notebook with critical variables defined at the beginning
• Apply zeros slip boundary conditions to classical smoothed solution

Things that we should figure out SHORT-TERM:

• Can we get displacements from relax at the elastic/VE interface and compare them with inferences?
• Can we generate a similar dip-slip case?
• What are the figures for the proposal?
  • 2-3 time steps?
  • Model geometry?
  • Predicted and recovered displacements
  • Estimated basal slip

Things we should look at LONG-TERM:

• Recovery with noise
• How far away from the fault should the mesh extend?
• Characteristic length scales for each eigenmode
• Make a more production notebook for
  • Autogenerate basal mesh from fault geometry?
  • Residual and basal slip roughness as a function of number of eigenmodes?
  • Remove explicit celeri dependence (We can just copy over the necessary functions over). This is so simple it would be nice to have it all in a single notebook.
  • Example with observational post-seismic from an earthquake or two (Hector, Tohoku)

[brendanjmeade]

[jploveless]

@jploveless please take a look at https://github.com/jploveless/relax_disloc/blob/main/relax_disloc_eigen_no_zero.ipynb.

Things that this does:

• Classical smoothed solution (without zero slip BCs)
• Eigenmode solution (with zero slip BCs)

This is really nice! I'm pretty fascinated by the differences between the estimated slip distributions, particularly the opposite polarity in the near-field tensile slip.

Things that I plan to do:

• Clean this up to make a nice demonstration notebook with critical variables defined at the beginning
• Apply zeros slip boundary conditions to classical smoothed solution

This sounds good. It'll be interesting to see if edge constraints on the smoothed solution have an impact on the overall shape/placement of lobes, and whether it brings the smoothed vs. eigen solutions any closer to each other.

Things that we should figure out SHORT-TERM:

• Can we get displacements from relax at the elastic/VE interface and compare them with inferences?

Yes — see note in the Google Drive email I just sent you.

• Can we generate a similar dip-slip case?

Currently saved on Shared Drive.

• What are the figures for the proposal?

  • 2-3 time steps?

I think we could get away with 2, one early and one late, focusing solely on the postseismic displacements rather than the total (coseismic plus postseismic).

• Model geometry?

I think the simple dipping is a good target here.

• Predicted and recovered displacements

I think the vector plots of horizontal surface displacements make most sense.

• Estimated basal slip

Yes. For compactness, I wonder if we could color triangles by slip magnitude, with vectors showing rake.

Things we should look at LONG-TERM:

• Recovery with noise

Definitely.

• How far away from the fault should the mesh extend?

This will be interesting to see, and also whether there's a difference between the yes/no edge constraints cases.

• Characteristic length scales for each eigenmode

And correlation of these length scales with layer thickness/viscosity parameters.

• Make a more production notebook for

  • Autogenerate basal mesh from fault geometry?

Seems straightforward enough. There are other files output by Relax that contain source fault depth information. That info is of course also contained in the input parameter file, which I could write a Python reader for in order to extract key parameters in a more automated way.

• Residual and basal slip roughness as a function of number of eigenmodes?

Yes, and I think it'll be interesting to see how this varies with VE parameters.

• Remove explicit celeri dependence (We can just copy over the necessary functions over). This is so simple it would be nice to have it all in a single notebook.

Seems straightforward. I was definitely being lazy in the initial notebook in grabbing the smoothing operator, etc. from celeri!

• Example with observational post-seismic from an earthquake or two (Hector, Tohoku)

Yes, and I think with these it'll be important to estimate slip both on a basal horizon and the source fault, so we can investigate true afterslip vs. VE that is mimicked by basal slip.

[brendanjmeade]

@jploveless Thanks for these comments! We're definitely on the same page. My two priorities are 1) a clean notebook and 2) A modification of the clean notebook for the dip slip case. Once we have the second, we can focus on figures for the proposal.

[jploveless]

Looking at the dip-slip case, it seems that the smoothing-based estimate yields slip distributions that are relatively compact, with positive and negative slip concentrations near the source fault. The eigenmode estimate (with 50 modes) shows something similar, but with high-magnitude, varying-sign "halos" closer to the edges of the basal mesh. Using fewer modes reduces these halos but results in a smoother slip estimate, with larger areas of lower magnitude slip but a generally poorer fit to the data.

I agree that some clean-up is important, but the notebook runs for the dip slip case just by changing the directory name!

[brendanjmeade]

Wow, that's awesome that we just need to point to the dip slip folder! So nice and it worked for me too. I'm seeing the same behavior that you are. In trying to investigate this, I'm trying to create a mesh with smaller spatial extent but a greater number of mesh elements. My thinking is that the higher wavenumber EVs may be mapping awkwardly to TDEs at low resolution. Setting the number of eigenvalues at 20 yields slip distributions that loom more similar.

However, I seem to keep crashing the kernel during the meshing stage when trying to specify mesh spacing < 0.5. For example, meshes[0].coords, meshes[0].verts = tri_horizon([-5, 5], [-5, 5], -2, 0.25) crashes the kernel with

Canceled future for execute_request message before replies were done
The Kernel crashed while executing code in the the current cell or a previous cell. Please review the code in the cell(s) to identify a possible cause of the failure. Click [here](https://aka.ms/vscodeJupyterKernelCrash) for more info. View Jupyter [log](command:jupyter.viewOutput) for further details.

Does this crash for you too? If so, thoughts on how to produce higher resolution meshes?

[jploveless]

Really interesting to see that a smaller mesh footprint seems to have more similar slip distributions across the smoothing and a smaller subset of eigenmodes! I think you might be right that the larger footprint (larger than the displacement grid) could be oscillating because of low resolution and/or low magnitudes of the peripheral part of the displacement grid.

About the crash, I think it's actually in cutde rather than in the meshing routine itself. I just commented out everything but the partials calculations from that cell, and it meshes with an element size of 0.25 no problem (and in 0.2 seconds) but I think it's a memory error in the partials calculation. Primarily, I think it has to do with the number of observation points in the full partials calculation (calculating disp_mat_full), because regular disp_mat, which uses every 10th observation coordinate, is fine (4.1 seconds). We don't actually need to do the full calculation; I was just doing it initially to exactly reproduce the .grd outputs of Relax.

[brendanjmeade]

Ahh that makes more sense. I'll note this as something to reorganize to create a clean starter notebook. Thanks!