Query about adding mechnochemical terms in thermal pressurisation model parameters set

[bowenyu1]

Dear QDYN developers, In our recent research, we found dynamic weakening could be greatly enhanced by mechanochemical reaction (e.g. flash heating promted by dehydration reaction and amorphization) for clay-bearing fault gouge. For further research, we want to see how this effect would contribute to the kinematics in a fault system. So I am writing to query if the mechanochemical terms (see the Eq. (6) and (7) and associated terms in the attached file) could be added in thermal pressurisation model parameters. Thanks. Best, Bowen Yu et al., 2023.pdf

[martijnende]

Hello and thanks for your suggestion. I had a quick look at the paper, which seems like a very complete body of work (congratulations!). It would be a rather big project to implement and test the equations that you refer to (and those that follow). I think a PhD student could spend 1-2 years on the implementation and the analysis alone. Depending on the degree of coupling between fluid pressure, temperature, and the chemical reaction kinetics, we could consider taking a few shortcuts. And if we can make it compatible with the current thermal pressurisation scheme (based on Noda's work), then the implementation would be a bit easier as the groundwork has already been done.

Do you have any plans for a project (with funding for a student or postdoc) to collaborate on this? There are a few other projects that are currently running in parallel that will bring big changes to the QDYN code base, which I will have to focus on first.

[bowenyu1]

Thanks for your reply. Considering the complexity as you mentioned, we decided to simply incorporate the flash heating model as well as the relationship between mechanical work and characteristic weakening temperature Tw (Eq. (1) in our paper （only 3-4 linear and quadratic equations）, and then see what the results look like. I am not sure if we could make this implementation through modifying the friction law in the script "friction.f90". Since the differential equations that describe the kinetics of chemical reactions have been abandoned, this may be simpler. I will appreciate it if you could give further suggestions to me. We have not made any plans for a project, simply for exploration at the current time. If we got one, we will contact you to collaborate on this in the future. We appreciate your interest. Bowen

[martijnende]

If you create a fork from this repository, then I can review some of the implementations you made to see if there's anything missing.

Note in relation to flash weakening: this formulation only works when the stress is computed from the velocity, and subsequenly the velocity is solved for with e.g. a Runge-Kutta scheme. This is currently the way QDYN works for rate-and-state friction, but in the future we will turn this around: the velocity is computed based on the current state of stress, and the stress is updated with the time-stepping routine. The reason for turning this around is the incorporation of creep with RSF (see #72)

Once we have turned the formulation around, the flash weakening approach no longer works, unless the rate-and-state parameters a and b are zero (so that the rate-and-state friction is constant). When b is zero, it is no longer possible to have unstable slip and so you will never reach the slip rate necessary for flash weakening.

So if you decide to use flash weakening, then be aware that your implementation will not work for future versions of QDYN.

[bowenyu1]

Thanks for your useful comments, I will try these implementations on my own first.

[martijnende]

Ok, I will close this ticket for now. Feel free to re-open and ask for advice if you need help. Good luck with the implementation!

[bowenyu1]

Dear Dr. Martijn van de Ende

These days I am attempting to incorporating a "new" friction law into QDYN that takes into account the regularized RSF law, flash heating, and the effects of mechanical work on weakening temperature Tw (see the equation set in the attached figure). The most difficult aspect of this approach was adding the differential equation "dmw/dt = tau * v" in order to calculate the real-time mechanical work mw and then Tw. In my modified version, the treatment of mw is identical to your treatment of the state variable theta. Could you tell me whether this implementation is correct?

I have successfully compiled the modified code, but an error titled occurs when I attempt to run the Single asperity (2D) tetorial code in python. The progression may halt at p.run(). I have tried to debug this in numerous ways, but have always failed. I would appreciate your assistance with that.

I have uploaded the updated filed assiciated with this issue. For your convenience, I have also categorized the positions I have made changes, as shown below:

1. In friction. f90: line18, 41-43, 78-79, 94-104, 113-114, 138-148
2. In problem_class. f90: line132 (neqs 3=>4), 117-118, 129-130, 151
3. In derivs_all. f90: line29, 38, 55, 99, 156
4. In solver. f90: line86, 161, 296
5. In utils. f90: line25, 29, 44, 59, 63, 76
6. In unittests_rsf. f90: line37, 95, 106
7. In unittests. f90: line171, 175. 188
8. In output. f90: line57, 93, 116, 372, 456, 458, 989
9. In input. f90: line58, 158
10. In pyqdyn. py: line114-119, 254, 432

Thanks for your time!

Bowen

modified src(bowen).zip

[martijnende]

Could you maybe create a fork of this repository and push your modifications to the forked repository? In that way I can more easily compare the changes that you've made to the code.

From a quick glance at friction.f90:79, you might be missing pb%theta_star in the equation. If I understood the equations correctly, your regularised RSF formulation should actually be the same as the one that was already implemented as case (2) in friction.f90, so I would just keep that as it is.

Then you also need to consider that in derivs_all.f90:152 you need to include the partial derivative of stress with respect to mw in the same way as for the state parameter. This line adds the radiation damping term to the stressing rate before being passed to the solver. If this radiation damping term is missing the mechanical work term (mw) then you will end up with a force imbalance that will crash the simulation.

I have successfully compiled the modified code, but an error titled occurs

Could you copy/paste the error message here? You can format it in a nice code style by using three backticks (3x ` ) to start a code block and three backticks to end the code block. Does the simulation run for a few steps or does it crash right away?

[bowenyu1]

Dear Martijn,

I have create a fork of your repository, the address is https://github.com/bowenyu1/qdyn.git.

The theta in my implementation refers to the fractional theta (i.e. theta/theta_star). I have tested these equations in a single matlab file, the output is fine.

For an easy start, the weakening temperature Tw is kept constant in the current forked repository. When all errors are overcome, I will introduce the evolution of Tw with mechanical work mw. I also need to think about how to construct the damping term in terms of mw. Could you give me any advice?

The codes were successfully compiled. But errors still exist when I attempt to run the Single asperity (2D) tetorial code in python. When I select the original RSF and aging law (both refer to case(0)), the error is titled:

Out[5]: -11

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:

Could not print backtrace: executable file is not an executable
#0  0x102c8e337
#1  0x102c8d313
#2  0x182419c43

When I select my newl RSF and state evolution law (both refer to case(4)), the error is titled:

Out[7]: -11

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:

Could not print backtrace: executable file is not an executable
#0  0x100e3a337
#1  0x100e39313
#2  0x182419c43

Thanks.

Bowen

[martijnende]

These segmentation errors imply that a variable that you're trying to access is not allocated. Can you successfully run the unit tests, or does it also crash when you try to run those?

As for the radiation damping terms: QDYN essentially solves the following equation

On the first line, the time derivative of the force balance ("stress equals strength") is written in terms of the partial derivatives of the friction coefficient µ. In your case, there is a partial derivative of µ with respect to Mw. The second line rearranges the various terms, putting the total derivative of the slip rate (v) on the left-hand side. What you will have to do is find an analytical expression of the partial derivative of µ with respect to Mw, and compute the rate of change of Mw (= tau * v).

[bowenyu1]

I can not run the unit test p.run(test=True, unit=True) either. It seems to crash without any steps. The error looks like the same :

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:

Could not print backtrace: executable file is not an executable
#0  0x102a2a337
#1  0x102a29313
#2  0x182419c43

As for the damping, I need to deduce it these days and then disscuss with you. Thanks.

[martijnende]

So it seems that somewhere during the initiation stage a variable has not been properly allocated. Segmentation faults are difficult to pin down, and in my experience the best way to debug this is by adding print statements at various places in the code, and see which ones get executed (and which don't). The default stdout channel is 6, so you can make a print statement like: write(6,*) "Test initialisation of solver".

Could you copy the complete output of p.run(test=True, unit=True), including all of the other messages? Perhaps that already tells us where the segmentation fault occurs. You should see at least something like Initiating QDYN unit test suite, which is being printed before the tests are being executed. If you don't see this, then the problem could lie in the initialisation of MPI.

[bowenyu1]

These are all the errors that I have seen. OK I will debug this following your suggestions. Thanks!

[bowenyu1]

Dear Martijn,

Recently I have tested a model with a depth profile which hold heterogeneous a-b, effective normal stress, pore pressure, effective normal stress and temperature (as shown in the figure below). Total modeling time is setted as 45 years. However, I could not get any results after an extremly long-time run (more than 24 hours), no matter the thermal pressurization is on or off. No erros occured. The situation did not change even though I tested the model on HPC. It seems like the calculation did not move forward after running at around a time of "0.4e-2 year".

Do you have any suggestions about that? Thanks.

Bowen

[martijnende]

I am currently working on a new release of QDYN (#84) that fixes lots of issues related to (effective) normal stress changes and MPI. Are you simulating a 3D normal fault by any chance? That is a known bug producing extremely small time steps and incorrect results.

[bowenyu1]

OK, I will try the new version of QDYN.

I am simulating a 1D fault in 2D medium. (MESHDIM=1).

[martijnende]

Eyup just mentioned to me that your normal stress is going to zero near the surface. This is problematic for the solver, as you are approaching a singularity (zero division). What happens if you offset your normal stress (something like $\sigma(z) = a z + \sigma_0$, with $z$ being the depth and $\sigma_0$ the effective normal stress at zero depth)?

[bowenyu1]

The code successfully run when I offset the normal stress using '''RNS_LAW = 0''' '''THETA_LAW = 1''' '''FEAT_TP = 0 or 1'''. However, the results look strange when I used '''RNS_LAW = 4''' '''THETA_LAW = 4''' '''FEAT_TP = 1''' (friction law considering flash heating, which was added by myself). It always stopped after a stress drop occured which is shown as the figure below: I have tried to added some limits on effective normal stress (1MPa when it is smaller than 1MPa) because of great weakening caused by flash heating. But it does not work. In fact, I have successfully tested and ran the code with '''RNS_LAW = 4''' '''THETA_LAW = 4''' under constant effective normal stress. Here is the results:

[martijnende]

There are probably additional dependencies on the normal stress that are not included in the partial derivatives (probably the dependency of the state parameter on the normal stress). See the comments starting here: https://github.com/ydluo/qdyn/blob/master/src/derivs_all.f90#L115

[martijnende]

Closing due to inactivity