HRF geometry files

[mallickrishg]

@brendanjmeade i looked at the new mht_geometry.csv files today, I think the format is perfect - 10 columns with a nice header makes everything easy to work with. Thank you so much for doing this!

Couple of points:

1. We need to y-coordinates to be flipped in sign. Currently I think you have the depth direction as $+y$ while the elevation direction is $-y$.
2. There is a column called 'name' in the csv file that I don't know what it is for. I was hoping you could clarify.
3. We need to discuss what the boundary conditions need to be for the various 'collection' items. Some are easy like 'topo' is traction 'bc' with bc_x=bc_y=0, similarly left and right edge are dirichlet bc's. But how we deal with 'base' and the 3 faults is not obvious, and we need to discuss this when you are free

[brendanjmeade]

Glad you can read this in.

1. We need to y-coordinates to be flipped in sign. Currently I think you have the depth direction as while the elevation direction is.

That's easy and I'll do that.

2. There is a column called 'name' in the csv file that I don't know what it is for. I was hoping you could clarify.

name doesn't have a deep meeting or purpose. Each collection is a conceptual grouping of elements "base", "mct", etc. name is just a local "index" into each element of collection. It doesn't have to be used but I thought might be useful in the sense that we might want to discuss something like the element named "n" in collection "topo". We could use the global index as well, but that would change for elements if the collections were specified in different orders. So I added this as. possible convenience.

3. We need to discuss what the boundary conditions need to be for the various 'collection' items. Some are easy like 'topo' is traction 'bc' with bc_x=bc_y=0, similarly left and right edge are dirichlet bc's. But how we deal with 'base' and the 3 faults is not obvious, and we need to discuss this when you are free

I totally agree and am similarly unsure. Some thoughts:

• base : Maybe eliminate this along with the left hand side vertical segment? We'd have to shorten the topo so that it terminates at the MFT.
• mbt: Some small slip rate (e.g., 5 mm/yr)
• ramp_mft: Spatially variable strike slip at element centroids. Are the rates from Judith's work or that by Sathiakumar? I'd guess that these are all 17 +/5 mm/yr.

• mct : Do not include for the first set of model runs. I recall that evidence for its Quaternary activity is equivocal.

  Once we agree on these, I'll implement them and send back a new .csv file.

[mallickrishg]

• base : Maybe eliminate this along with the left hand side vertical segment? We'd have to shorten the topo so that it terminates at the MFT.

I don't think we can truncate topography. Ideally we need to extend topography arbitrarily far away from the faults. However, I think we can get ride of base and left_edge. Regarding mct - I agree we can set it to have 0 slip, and so it acts only as a passive observation point for stress changes.

For mft and mbt we need to divide up the ~15mm/yr of convergence (which might become something like 16-18mm/yr of slip rate for those fault dips) into the 2 faults somehow. We could actually just do a bem solve assuming those faults are traction-free surfaces accommodating the overall convergence boundary conditions applied at the right_edge.

Additionally we need to come up with how we label bc_type. I think there could be 3 types

• u - displacement bc (will have associated displacement kernels)
• t - traction bc (will have associated traction kernels)
• s - imposed slip bc (this is where the linear operator from get_matrices_slip_slip_gradient comes in)

[mallickrishg]

@brendanjmeade I think bemcs is nearly ready to roll and do some real HRF problems. In terms of code we still have a couple of things we need:

• we need to allow input files to have mixed slip and traction boundary conditions for a given fault patch
• need to specify if tractions or slip are in global (x,y) or local (s,n) coordinates for a given patch

I look forward to hearing your thoughts on this.

Most important thing I found was the condition number of the linear operator goes from $6\times 10^6$ when we have traction boundary conditions in a closed polygon down to $10^3$ when we extend topo by just 1 element outside the closed polygon!

Attached is a working example showing the displacement field for the most recent file you shared (with an additional element extending topo to $x>0$. Currently it's in a branch called 'hrf-bem'. I will merge it later - still testing some things. https://github.com/brendanjmeade/bemcs/blob/hrf-bem/examples/hrf/example_hrf_bem.ipynb

[brendanjmeade]

@mallickrishg Wow! I can believe and I'm so happy to see it.

• Regarding condition numbers. We can always precondition the linear operator to get a better condition number. I shared some basic examples in: https://github.com/brendanjmeade/bemcs/blob/main/notebooks/example_matrix_preconditioners.ipynb.

  We’ll probably need to get used to this once we approach larger solves where we’ll do interactive solutions rather than direct ones. We can pull that and the h-matrix solver right out of the celeri block modeling code.

• Regarding boundary conditions. How about this: six types of possible boundary conditions:

  • u_local
  • u_global
  • t_local
  • t_global
  • s_local
  • s_global

  We can specify in a .csv geometry/BC file and then the BEM code can process as appropriate.

  As to the second part of the problem, "mixed boundary conditions for a given fault element”. As I understand it, this type of thing is complicated because it requires a relative weighting of the two types of boundary conditions. There is an equation for this in the definition of a Robin boundary condition that augments the linear system, but it introduces a new tunable parameter. I’ve never studied these before I’ve only ever discussed them. Omar Ghattas has done these for glacier problems, and he’s got some optimization technique for this. I’m happy to consider this problem but think we should focus on classical mixed BC problem first and leave the Robin problem till we’ve really got the fundamentals down. What do you think?

[mallickrishg]

I like your suggestion, and I agree we should keep the bc labels in the input file itself.

Regarding the mixed conditions, what I meant is not Robin bcs, although that's what we need for a proper friction problem. But that's for a later day. I was instead talking about allowing the x and y bcs to not have to be the same type. Basically that would mean 2 columns for bc_type and 2 columns for bc_x,bc_y.

For example for ramp_mft type of elements we specify each of the elements to have

• 's' type bc along the fault-normal direction with value 0
• 't' type bc in shear direction with value 0

[mallickrishg]

I was thinking about the preconditioners problem yesterday.

I don't know how a preconditioner can help solve the problem when I think the issue is that closed polygon geometries can be translated freely while still satisfying traction boundary conditions, and to a limited extent displacement boundary conditions.

Column or row preconditioning might reduce the condition number, and we may need it anyway when we are dealing with really large matrices, but that may not solve the more fundamental issue of unestimated translational parameters. By extending the mesh outside of a closed polygon to create 'open' boundaries (which are forced to have 0 slip), we anchor the solution. For the hrf example, by extending topo on both sides we solve this problem!

[brendanjmeade]

I don't know how a preconditioner can help solve the problem when I think the issue is that closed polygon geometries can be translated freely while still satisfying traction boundary conditions, and to a limited extent displacement boundary conditions.

I understand now. Maybe the Crouch and Starfield book have suggestions in their fictitious force section where they only use traction BCs?

[mallickrishg]

This is a sort of gap in my understanding of BEM - I don't actually know when is a good idea to use displacement discontinuity vs fictitious force method especially when we only have displacement and traction boundary conditions.

If you gave me slip BCs then I know that we have to use displacement discontinuity and I don't think it is even possible to solve the same problem with fictitious forces. However, if we only have forces or loads or equivalent tractions as BCs, then I think it may be possible to solve the same boundary value problem with either method (albeit with different coefficients). I was hoping you have some insight into this?

[brendanjmeade]

Yeah traction only BC problems can translate arbitrarily. I try to stay away from thinking about them because of this issue.

[mallickrishg]

Btw - I have creating a file with all the automatic labelling functions in it : https://github.com/brendanjmeade/bemcs/blob/hrf-bem/examples/labeling.py

I was wondering if I should add this to the bemcs monolith? What do you think?

[brendanjmeade]

This is awesome and I definitively think it should be added to the monolith and exported in__init__.py. Thanks for the thumbs up on the boundary condition idea too. I'll get that implemented in the input .csv file. I really like your suggestion to have a lot of BC options in the input file.

[mallickrishg]

When you have time, could you send me the new mht.csv file with the new boundary conditions and topography extended in both directions? I will setup the notebook to deal with the new input file (hopefully this is the format we will stick with), and test and push the labeling/node handling functions to the bemcs monolith.

[brendanjmeade]

Absolutely! It's on my schedule for Wednesday morning. Pardon any delay.

[mallickrishg]

We have our first proper Himalaya example ready to roll! I have merged my working branch into 'main', and so you can play around with the notebook: https://github.com/brendanjmeade/bemcs/blob/main/examples/hrf/example_hrf_bem.ipynb

The really nice part about implementing 'local mixed BCs' is that we have purely fault-parallel slip at the central nodes of each element, even if the overall geometry is non-planar - which is what all geologists seem to enjoy talking about! The additional utility of these kinds of boundary conditions is we can actually attempt static-friction boundary conditions (if we wanted) i.e.,

1. BC1: $\tau_s - f\taun = 0$ which we can rewrite as a linear equation using both traction kernels as: $\left[K{\taus} - fK{\tau_n}\right].[s_s,s_n]^T$ where $K$ is a traction kernel
2. BC2: $s_n = 0$

not something we need to implement or play with right now, but is the 'correct' way to address problems about non-planar or rough faults. I found some interesting ideas in the following papers:

1. https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2000JB900241
2. https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2011JB008624
3. https://www.sciencedirect.com/science/article/pii/S0191814118301597

[mallickrishg]

For next steps re the paper, I was thinking of 3 input files.

1. The current mht geometry with boundary conditions
2. The brittle locked interface model, with ductile part shear traction free
3. Part of the interface is locked and we impose the 2015 Gorkha event

[brendanjmeade]

We have our first proper Himalaya example ready to roll! I have merged my working branch into 'main', and so you can play around with the notebook: main/examples/hrf/example_hrf_bem.ipynb

This is fantastic! On Friday I'm going to carve out time to plot up the internal stresses and strain energy.

The really nice part about implementing 'local mixed BCs' is that we have purely fault-parallel slip at the central nodes of each element, even if the overall geometry is non-planar - which is what all geologists seem to enjoy talking about!

NGL this is pretty brilliant and I totally wouldn't have thought of this. This is the best project ever.

The Ritz and Pollard paper and the Davis papers are new to me. I've love to revisit these. They all lean on the old Saucier 1992 paper which showed that even small variations in geometry can have big effects on off fault stress orientations.

[brendanjmeade]

2. The brittle locked interface model, with ductile part shear traction free

What part do we want to define as ductile? Should we consider extending the ramp_mft structure farther to the north?

3. Part of the interface is locked and we impose the 2015 Gorkha event

For the Gorkha event, I imagined that we'd do it with just the Gorkha slipping faults and then the free surface. Nothing else would even exist. Maybe that's the same as saying the interface is locked?

[mallickrishg]

Regarding locked and ductile/creeping, we can just use any of the interseismic coupling models for this : they all show the locked section that is atleast 100km inland from the fault trace. We can those parts of the fault (s_s,s_n = 0,0). The rest of the mht is traction free, s_n = 0.

For the Gorkha event, I imagined that we'd do it with just the Gorkha slipping faults and then the free surface. Nothing else would even exist. Maybe that's the same as saying the interface is locked?

Yup, they are equivalent. We will set slip to be 0 everywhere except on the Gorkha rupture and topo, thereby making it so that those faults don't exist.

[mallickrishg]

@brendanjmeade - I wanted to update you that I have made some changes to the bemcs monolith:

1. I realized that in some cases the reordering of elements in standardize_els_geometry() makes it impossible to actually run the code, especially when the domain is a box! I have now added a flag called 'reorder=True/False' to give the user some flexibility (the node labelling scheme will still throw an error if the circulation pattern is incompatible). The best example of the need for this is in the 'fault in a box' problem - bonus: here I show how one can handle static friction as a boundary condition - yaaay!
2. Turns out there were bugs in dealing with node labelling of closed objects - they are now fixed!

Do you want me to go ahead with creating the 3 Gorkha experiments? I don't want to go too far ahead because I want to give you time to sanity check some of the more recent features I wrote.

[brendanjmeade]

@mallickrishg

• Regarding (1.): that's interesting probably something that I should have anticipated. Apologies if the order scheme I proposed caused too many troubles over the last month.

• Regarding (2.): Nice and finding these things is a part of the process!

• I think moving ahead with the Gorkha experiments is a good idea. I've finished some editing of the theory part of the text now and feel good that we agree on it1

[mallickrishg]

@brendanjmeade A notebook with the following calculations for the HRF example is ready here: https://github.com/brendanjmeade/bemcs/blob/main/examples/hrf/experiment_hrf_eqcycle.ipynb

It contains calculations for the following scenarios: (1) steady state/long-term creep on the plate boundary (2) late interseismic shallow locking (3) Gorkha partial rupture coseismic

[mallickrishg]

@brendanjmeade constant slip BEM also done: https://github.com/brendanjmeade/bemcs/blob/main/examples/hrf/example_hrf_constantslip_vs_bemcs.ipynb

I've shown the calculations for the Gorkha event case, you can swap it out for any other case by simply changing the name of the file and it should work. The solution, in terms of the recovered coefficients, are mostly similar between the "smooth" formulation and the constant slip approach. This is to be expected i guess. The big difference is once you plot up the internal stresses!

[brendanjmeade]

Awesome and thank you. The notes are useful to me too! My goal is to get this work for the interseismic case and show the internal strain energy for both cases. That should be pretty dramatic