Slip azimuth constraints

[brendanjmeade]

Rishav Mallick (@mallickrishg) mentioned that he might be interested in slip azimuth constraints. I've never used these before. The basic idea is to provide a constraint on the relative magnitudes of ($t\mathrm{e}$, $t\mathrm{n}$) on some TDE element or ($s\mathrm{ss}$, $s\mathrm{ts}$) on some fault segment. The elements:

Formulation

TDE azimuth constraint. We probably want to specify the angle, $\alpha$, at some location. An obvious way to specify this is as,

$\tan \alpha = t\mathrm{n} / t\mathrm{e}$.

This is non-linear in the estimation parameters, $\mathbf{t}$. So we need to re-write this as something linear. How about:

$0 = t\mathrm{n} - t\mathrm{e} \tan \alpha$.

which is now linear in $\mathbf{t}$. We can just add this as an extra row of constraints to operator in celeri.

Slip rate azimuth constraint. Something similar to above but it needs some rotations to get into east and north coordinates from strike-slip and tensile_slip coordinates. Not sure.

File format

TDE azimuth constraint. A .csv file typically named something like slip_azimuth_tde.csv.

mesh_name, tde_index, azimuth, azimuth_uncertainty

TDE azimuth constraint. A .csv file typically named something like slip_azimuth_segment.csv.

segment_name, azimuth, azimuth_uncertainty

[mallickrishg]

@brendanjmeade I think what you have written out in the linear problem is correct, but there might be an easier way to write this. The dataset you get from focal mechanisms are the strike ($\phi$) and rake ($\lambda$) angles, which one can convert to slip vector azimuths ($\phi - \lambda$) OR convert to a new rake angle ($\lambda'$) given the strike ($\phi'$) of the block boundary fault of interest. If the focal mechanism and block boundary strikes match exactly, then one can use the focal mechanism rake directly.

So, for the 'observed' rake angle for any fault patch (rectangle or TDE), the only equation you need is the slip rate in the rake perpendicular direction is 0 i.e., $[R(\lambda' \pm \pi/2)] [\omega] = 0$ where $R$ is the relevant rotation matrix.Alternatively, if you already have equations that give you slip rates in the strike-slip and dip-slip directions then you just need to rotate those basis vectors to rake-parallel and rake-perpendicular and zero out the latter.

[brendanjmeade]

@jploveless Thoughts on this? I have no immediate plans to work on this but want to document ideas about how to best do this.

[jploveless]

I really like the approach @mallickrishg suggests, of rotating the strike-slip and dip-slip partials to the rake-parallel and rake-perpendicular basis vectors by way of the rake itself. I think (but haven't yet drawn it out to clarify to myself!) that this could eliminate a potential $\pm \pi$ ambiguity that could exist in a direct approach of setting $0 = t_n - t_e \tan α$.

A related constraint on TDE slip is one that I had implemented in Blocks: setting coupling equal to 1. I feel this is a good approach at the lateral edges of TDE meshes, because the adjacent rectangular segments are by definition fully coupled. The constraint on these edge-lining TDEs therefore sets the estimated slip rates to be equal to projected relative block motion rates, which uses a function that transforms relative block motion onto the TDE geometry. I mention this here because I could see the potential for a TDE constraint in celeri being "constrain TDE slip rake to be consistent with relative block motion" (without necessarily enforcing full coupling). I don't think this is always a good idea (because spatially variable slip rates should require spatially variable rakes) but I could see it being a parameter worth testing in some cases.

[brendanjmeade]

Thanks @jploveless your coupling approach in mind to you think we need another way to specify an azimuth prior on TDEs or is what we have sufficient. Either way, what are your thoughts on where to store this information? As additional fields in a *mesh_parameters.json file? Similarly we could store fault segment slip azimuth priors in a *segment.csv file. That would require updating celeri_ui but I'm willing to do that. Thoughts welcome!

[jploveless]

Sorry — I completely spaced on that part of the thread.

I think this is worth thinking about, because it could prompt a bigger restructuring of how a priori fault information is stored. We currently specify a priori slip rates for block boundaries as part of the *segment.csv file. I see this as sensible because it keeps the specified slip rate as a property of the segment, and therefore it would make sense to specify rakes/azimuths in the same way. On the other hand, I don't think that specifying a priori slip rates/rakes/azimuths for TDEs can be cleanly done in a *mesh_parameters.json file, because those contain properties of the entire mesh, not individual elements. In Blocks, specifying a coupling fraction as a priori information on TDEs is a little clunky, requiring a simple, separate text file that specifies element index and slip rate or coupling fraction, but it works, and we could easily implement something similar in celeri.

One new approach could be to move all a priori slip information to a new *.csv file whose lines contain geometry/location and slip information. I see a priori slip constraints as an additional dataset, like geodetic data. I could imagine slip information being specified just in terms of latitude/longitude (plus rate, rake, or azimuth, with some flag indicating the type of constraint). A helper function could then map the slip information onto the correct segment or TDE, by finding the segment/TDE centroid that is closest to the constraint location. I could see this being a more natural way of using focal mechanism azimuth/rake constraints, although I think it loses some of the convenience of celeri_ui of specifying segment properties.

What do you think?

[brendanjmeade]

I think this is worth thinking about, because it could prompt a bigger restructuring of how a priori fault information is stored. We currently specify a priori slip rates for block boundaries as part of the *segment.csv file. I see this as sensible because it keeps the specified slip rate as a property of the segment, and therefore it would make sense to specify rakes/azimuths in the same way.

Agreed that this isn't bad. All fault segment information in one place.

On the other hand, I don't think that specifying a priori slip rates/rakes/azimuths for TDEs can be cleanly done in a *mesh_parameters.json file, because those contain properties of the entire mesh, not individual elements.

We could add fields to *mesh_parameters.json that contain information about things on specific TDE indices. This has usability and elegance level of -1000 but we could do it.

In Blocks, specifying a coupling fraction as a priori information on TDEs is a little clunky, requiring a simple, separate text file that specifies element index and slip rate or coupling fraction, but it works, and we could easily implement something similar in celeri.

Sensible and easy...I like it. Something like a .csv file with:

mesh_name, mesh_element_index, ss_slip_rate, ss_slip_rate_sigma, ss_slip_rate_flag, etc. 

One new approach could be to move all a priori slip information to a new *.csv file whose lines contain geometry/location and slip information.

First of all this is frightening. Second of all it's definitely worth discussing.

I see a priori slip constraints as an additional dataset, like geodetic data. I could imagine slip information being specified just in terms of latitude/longitude (plus rate, rake, or azimuth, with some flag indicating the type of constraint). A helper function could then map the slip information onto the correct segment or TDE, by finding the segment/TDE centroid that is closest to the constraint location. I could see this being a more natural way of using focal mechanism azimuth/rake constraints, although I think it loses some of the convenience of celeri_ui of specifying segment properties.

I definitely can see that, but it introduces a layer of ambiguity that doesn't appeal to me. Specifically the closeness idea. Currently, we specify this exactly and taking a step away from saying exactly where is a big step for uncertain gain. Certainly willing will ponder and happy to change my mind.

One other possibility for the TDE case (which is likely too radical) is to move away from the .msh format and specify meshes in a .csv file. Yep crazy...but just for the sake of argument. Consider a format like this:

mesh_id, x1, y1, z1, x2, y2, z2, x3, y3, z3, ss_rate, ss_sigma, ss_flag, etc.

In this case, all mesh geometry, slip rates, and prior constraints could all be stored in a single hand-editable file for each mesh. If one wanted to, multiple meshes could even be concatenated into the same file. There's lot to like about that. On the downside, this file would have to be created from a .msh file or some other meshed file format...so maybe it's just passing the buck? Not sure, but I thought I'd bring it up. Thoughts (including vehement rejections) are welcome.

The easiest thing to do is probably just to create just add slip azimuth constraints to *segment.csv files and add a new file for TDE constraints...but it's also fun to think about the right thing to do.

[jploveless]

I think this is a really clear layout of the pros and cons of these approaches. I agree that the "assign slip constraints by proximity" could be fraught with ambiguity, and I do think it's important to make sure the constraints are being applied to specified segments/elements as current practice does.

A twist on your suggestion to move to .csv for TDEs, which I see as essentially an equivalent to a *segment.csv file, is to dive deeper into the capabilities of .msh files. There are some optional fields within Gmsh files that I think could be used to assign an arbitrary number of properties, which really stem from its usage as a pre-processor for using the generated meshes in a FEM. I'm not sure how these optional fields play with the existing Python IO tools, but it might be worth exploring. It's also worth noting that a single .msh file can contain multiple distinct entities, and so it is also possible to use a single file to hold information about multiple meshes, rather than using multiple files.

All this said, I agree that the easiest solution is your final suggestion. Part of this could be to write a writer that parses geometry information from .msh files and constraint information from a separate .csv file and writes a combined .csv file using the format you suggest.

[brendanjmeade]

A twist on your suggestion to move to .csv for TDEs, which I see as essentially an equivalent to a *segment.csv file, is to dive deeper into the capabilities of .msh files. There are some optional fields within Gmsh files that I think could be used to assign an arbitrary number of properties, which really stem from its usage as a pre-processor for using the generated meshes in a FEM. I'm not sure how these optional fields play with the existing Python IO tools, but it might be worth exploring.

That's a great idea and probably worth exploring. I have sort of lingering worry about support for reading .msh files in python because the developer of the meshio libraries has decreased their focus on that work. There are a large number of outstanding issues, the repo hasn't been updated in a year, and the developer isn't active on Discord anymore. Having said that we can probably roll our own if the need arises.

All this said, I agree that the easiest solution is your final suggestion. Part of this could be to write a writer that parses geometry information from .msh files and constraint information from a separate .csv file and writes a combined .csv file using the format you suggest.

Slip azimuth as a constraint in a segment file seems reasonable and is low-cost to implement.

Slip azimuth as a constraint on TDEs is still open for discussion.

[brendanjmeade]

Returning after a long hiatus. Here's my current input file proposal for applying slip rate constraints to meshes. It simply adds fields to a *mesh_parameters.json file. The files are:

1. indices of elements with strike-slip constraints
2. strike-slip constraint values
3. strike-slip constraint uncertainties
4. strike-slip constraint weights (not sure if this is useful)
5. indices of elements with dip-slip constraints
6. dip-slip constraint values
7. dip-slip constraint uncertainties

8. dip-slip constraint weights (not sure if this is useful)

   and each of the fields is a list. For a concrete example, see the eight lines at the bottom of the json below (file attached as well). This approach is not convenient for large numbers of constraints, but it is very straight forward. I'm approaching an application of this, so I'd welcome feedback from @jploveless, @mallickrishg, @Emilycm-97

[
  {
    "mesh_filename": "../data/mesh/cascadia.msh",
    "smoothing_weight": 1e18,
    "n_eigen": 20,
    "top_slip_rate_constraint": 0,
    "bot_slip_rate_constraint": 1,
    "side_slip_rate_constraint": 1,
    "top_slip_rate_weight": 1,
    "bot_slip_rate_weight": 1,
    "side_slip_rate_weight": 1,
    "a_priori_slip_filename": "",
    "ss_slip_constraint_idx": [2, 43],
    "ss_slip_constraint_rate": [0, 0],
    "ss_slip_constraint_sig": [1, 1],
    "ss_slip_constraint_weight": [1e3, 1e4],
    "ds_slip_constraint_idx": [2, 43],
    "ds_slip_constraint_rate": [5, 5],
    "ds_slip_constraint_sig": [1, 1],
    "ds_slip_constraint_weight": [1e2, 1e2]
  }
]

and this can be read with:

import json
import rich
mesh_params = json.load(open("cascadia_mesh_parameters_with_constraints.json"))
rich.print(mesh_params)

to produce a dictionary:

[
    {
        'mesh_filename': '../data/mesh/cascadia.msh',
        'smoothing_weight': 1e+18,
        'n_eigen': 20,
        'top_slip_rate_constraint': 0,
        'bot_slip_rate_constraint': 1,
        'side_slip_rate_constraint': 1,
        'top_slip_rate_weight': 1,
        'bot_slip_rate_weight': 1,
        'side_slip_rate_weight': 1,
        'a_priori_slip_filename': '',
        'ss_slip_constraint_idx': [2, 43],
        'ss_slip_constraint_rate': [0, 0],
        'ss_slip_constraint_sig': [1, 1],
        'ss_slip_constraint_weight': [1000.0, 10000.0],
        'ds_slip_constraint_idx': [2, 43],
        'ds_slip_constraint_rate': [5, 5],
        'ds_slip_constraint_sig': [1, 1],
        'ds_slip_constraint_weight': [100.0, 100.0]
    }
]

cascadia_mesh_parameters_with_constraints.json

[jploveless]

I think this is a nice solution. Specifying these constraints in the mesh parameter file feels natural and organized. I agree that it could get cumbersome to apply too many constraints, but I don't feel like that will be a common exercise.

One potential shortcut could be to allow items like ss_slip_constraint_rate to either have the same number of elements as the corresponding *_constraint_idx field, or be a scalar, in which case it's assumed that the same rate, sigma, and/or weight applies to each constraint. But that would likely only be useful for larger numbers of constraints than the 2 you're showing in the example here.

I look forward to seeing your application of this!

[brendanjmeade]

@jploveless

One potential shortcut could be to allow items like ss_slip_constraint_rate to either have the same number of elements as the corresponding *_constraint_idx field, or be a scalar, in which case it's assumed that the same rate, sigma, and/or weight applies to each constraint. But that would likely only be useful for larger numbers of constraints than the 2 you're showing in the example here.

Thanks for reviewing this. It'll be good to have a decent mechanism to get these data into celeri.

I like your shortcut idea too. Implicitly broadcasting a single value to all named indices is pretty clever!

[jploveless]

979db61 (or rather ba5ca3f, which includes a notebook I had forgotten to add to the previous commit) starts to chip away at ideas mentioned in this issue, but not yet to the point of allowing azimuth constraints. But, specifying slip rate and coupling constraints is now possible, within the mesh_param.json file using @brendanjmeade's suggested syntax. A new notebook doesn't include much new stuff but shows a representation of the projected block motions as the last code cell.

• Mesh edges (top, bottom, side) can have no constraints, a constraint of full creep (forced to be zero slip deficit rate), or full coupling (where slip deficit rate in both directions is forced to be the same as the projected relative block motion)
• Arbitrary elements can have slip rates assigned in either/both directions
• Arbitrary elements can have any coupling fraction assigned. This forces the estimated slip deficit rate components to be that fraction times the projected relative block motion rate components.

As of now, individual weighting and sigmas are not honored. These can be specified in the mesh_param.json file but currently the weighting of triangle element slip constraints is controlled by tri_con_weight. I'll try to integrate those individual weight/sigmas soon, as well as work on @mallickrishg's suggestion for slip azimuth constraints.