Sir, I am using the subduction3d example in PyLith 2.2.2 to perform elastic and viscoelastic simulations in step06. Then, I use step07a and 07b to compute the Green’s functions. However, I encountered a problem: although the surface GPS displacements resulting from different models in step06 are clearly different, the inversion results using the Green’s functions from step07a and 07b turn out to be the same. Could you please help me identify what the issue might be?
file.zip (402.6 KB)
I set the elastic and viscoelastic parameters in step06 to obtain the forward modeling results under both elastic and viscoelastic conditions in order to compute the Green’s functions. Is this approach incorrect? I’m truly sorry to bother you again, but due to the urgency of the situation, I’m reaching out once more. Thank you sincerely for taking the time to read this!
I do not understand what you are trying to do. Please provide a more complete description of what you are trying to do, using diagrams and equations as necessary.
What do you mean by “although the surface GPS displacements resulting from different models in step06 are clearly different, the inversion results using the Green’s functions from step07a and 07b turn out to be the same”? What are the “different models in Step6”? Be more specific about exactly what inversion you are doing. What are the inputs? What are the outputs? What do you mean by “the same”?
Sir, this is my thought: in your example, you showed that we can calculate surface GPS displacements based on prescribed fault slip using step06, and then generate the cgps_synthetic_displacement.txt file using make_synthetic_gpsdisp.py. After that, we compute the Green’s functions via step07a and step07b, which are ultimately used for inverting fault slip rates. In my case, the inversion data is lindsey4.txt, and the inversion script is 1slip_invert.py.
Previously, I have only tested the elastic model (in step06). Now I want to invert fault slip using a viscoelastic model, in order to consider the influence of viscoelastic relaxation, and then compare it with the elastic results (using the same inversion data and script as before). I found that the cgps_synthetic_displacement.txt files generated by make_synthetic_gpsdisp.py for the elastic and viscoelastic models are significantly different. This indicates that the outputs from step06 are clearly different (note: the elastic model uses a time period of 1 year, while the viscoelastic model uses 10,000 years; I compared the results based on scaling in the displacement file, and also compared between two viscoelastic models directly).
Given this, I believe the Green’s functions computed in step07a and step07b should also differ between elastic and viscoelastic models. However, after running the inversion with lindsey4.txt and 1slip_invert.py, I found that the resulting fault slip distributions for all three models are the same. This seems incorrect to me. Could there be a problem with my procedure?
By the way, although I used mat_elastic.cfg in my viscoelastic model, the material definitions inside are actually viscoelastic.
When you invert the displacements for the simulation with viscoelastic deformation, what time snapshot from the output are you using? If you use the first time step in the simulation output for the displacements in the inversion, then it likely matches the elastic deformation, so the results should be the same as in the elastic deformation. If you use the final (last) time step of the displacement output, then the results should be different. You should be able to debug the inversion script to verify which time step in the output you are using for the displacements that you are inverting.
As you can see from the file above, I used the 10,000-year timescale throughout the process in make_synthetic_gpsdisp.py. Isn’t the generated data intended for computing Green’s functions? Do you mean that I should make the corresponding time adjustment in 1slip_invert.py?
in step06:[pylithapp.problem.formulation.time_step]
Define the total time for the simulation and the time step size.
total_time = 10000.0year
dt = 0.8year,and in make_synthetic_gpsdisp.cfg:# We select the last time step for this problem.
time_step = 12500
Could you kindly help me check where there might be a problem in my file? I’d be very grateful!
Dear baagaard, Could it be that I didn’t set the time interval for computation in step07a and step07b, which caused all the Green’s functions to correspond to only the first timestep? From my understanding, the cgps_synthetic_displacement.txt obtained in step06 is used to compute the Green’s functions, so I thought my previous approach was actually calculating the Green’s functions for the last timestep. Is it that I’m supposed to set the timestep in step07a.cfg and step07b.cfg instead?
I am still having difficulty understanding what you are trying to do. The PyLith Green’s function problem type is designed for static Green’s functions. So it will give you the static impulse response at all station locations resulting from slip impulses on the fault. It is not designed for viscoelastic materials (which would involve computing the quasistatic viscoelastic response for each slip impulse, not the static elastic response).
If you have actual observations you want to invert for static fault slip, then you do not need “Step 6”, which computes fake observations. Instead, you feed your actual observations to the inversion script, and it will invert for the static slip.
If you are trying to compute viscoelastic Green’s functions, then you will need to setup your own workflow for the inversion. You would compute the Green’s functions by running a separate forward PyLith time-dependent quasistatic simulation for each slip impulse that you create on the fault. You will also have to create your own inversion script, because the one provided in the subduction-3d example is designed for static slip.
I understand what you mean. I would like to ask whether there is a specific template for “a separate forward PyLith time-dependent quasistatic simulation for each slip impulse.” In my step07, under:# ----------------------------------------------------------------------
faults
----------------------------------------------------------------------
[pylithapp.problem]
interfaces = [fault_sg]
To create slip impulses for Green’s functions we use a special fault type.
interfaces.fault_sg = pylith.faults.FaultCohesiveImpulses
[pylithapp.problem.interfaces.fault_sg]
id=10
Impulses for reverse slip.
Note that it is possible to apply both left-lateral and updip slip
(impulse_dof = [0,1]), but we separate the impulses into two problems.
impulse_dof = [1]
how should I modify this?
If I only need to calculate the Green’s function after reaching the steady state, does it mean I just need to modify the material parameters in the original elasticity calculation and add a total time simulation?
Dear Sir,
I’m very sorry to bother you again. I’ve tried many approaches and also reviewed the PyLith user manual. However, it only provides guidance on computing static Green’s functions. As you mentioned, computing viscoelastic Green’s functions involves time-evolving fault slip impulses, but I’m still quite confused about how to correctly configure the .cfg files for such calculations.
I noticed your suggestion to use TimeDependent instead of greensfns, and that the spatial database should vary with time. Based on this, I attempted to set up the necessary files as shown below, but the simulation still failed. Could I kindly ask you to take another look to help identify where the issue might be?
Right now, my main confusions can be summarized into the following three points:
- How to assign different impulses at different times on each point of the fault surface — that is, how should the spatial database be constructed in this case?
- How to correctly configure the
.cfgfiles — similar tostep07a.cfgandgreensfns.cfg— to ensure the time-dependent simulation runs successfully? - I aim to compute Green’s functions for a viscoelastic model after reaching steady state. In this case, would it be acceptable to compute Green’s functions only at the final time step, and then use your
slip_invert.pyscript for inversion?
Thank you again for your patience and guidance.
greensfile.zip (8.9 KB)
If you want to compute viscoelastic Green’s functions, then you cannot use the GreenFns problem or the FaultCohesiveImpulses. You need to follow a normal viscoelastic forward simulation like examples/subduction-3d/step02_cosesimic.cfg but replace the slip distribution with a slip that is an impulse. You will have to create a spatial database that corresponds to a small slip patch of slip surrounded by zero slip, and do this for each impulse you want in your inversion. If you want the displacements of the impulses to correspond to the final completely relaxed state, then run each simulation to a total time that is several Maxwell times. In the inversion, you would use the final displacements.
Sir, I still don’t quite understand how to implement “You will have to create a spatial database that corresponds to a small slip patch of slip surrounded by zero slip, and do this for each impulse you want in your inversion.” Is it similar to what is described in hex8/step021.cfg?
I’m not sure if I understand this correctly: I need to prepare many spatial files, each of which defines slip in only one specific region while the rest have zero slip. Then, I need to run the simulation once for each of these spatial files to generate the corresponding number of output files. In the inversion process, I will read the data from all these output files to construct the Green’s function matrix for inversion.
Dear baagaard
I tried it this afternoon, and please find my attachment below. My idea is to divide the fault into a 6×150 grid. In the first simulation, I assign a left-lateral strike-slip value of 1 to a single grid cell, with all others set to 0. I then repeat this process for each cell to obtain output1 through output900, recording the fault slip and the displacement at GPS stations over a 100-year time scale. By reading the slip distribution and the corresponding GPS displacement after 100 years for each cell, I aim to construct our Green’s functions.
However, I’m not sure whether this approach is correct, and I’m also uncertain if the files in the attachment are accurate. Would you mind taking some time to review them? I would greatly appreciate it.
I also have a few questions:
- I have only calculated strike-slip components. Should I also simulate 900 thrust-slip components?
- Is there a simpler way to achieve this without running the simulation 900 times and preparing so many spatial database files? For example, is it possible to run the model once and read from a single spatial file instead?
- How should I modify
slipinvert.pybased on this approach? - In the above simulations, I used the
simplGridDBspatial database, because usingsimpleDBcaused an error. However, I’m not sure what the problem was. Will this affect the accuracy of the Green’s function calculation? - In my simulations, I noticed that multiple files are generated, such as
step07a-leftlowercrust.h5_t2997972000000002.vtk, but they only contain information for theleftlowercrustmaterial. I’m not sure why the other materials are not included.
Thank you very much for your time and kind attention. I wish you all the best!
viscoelastic.zip (3.0 MB)
I have only calculated strike-slip components. Should I also simulate 900 thrust-slip components?
If you want to allow spatial variation in the rake angle, you need both components, just like step07a and step07b in examples/subduction-3d.
I strongly recommend considering the size of the fault patches. Do the observations constrain slip at 900 points? There are studies that demonstrate reducing the resolution as a function of depth, for example https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2004GC000841.
Is there a simpler way to achieve this without running the simulation 900 times and preparing so many spatial database files? For example, is it possible to run the model once and read from a single spatial file instead?
There is not a simpler way with the current PyLith application. It is possible to create a high-level Python script to call the PyLith application.
How should I modify
slipinvert.pybased on this approach?
You will have to loop over all of the simulations producing Green’s functions to assemble the matrix with Green’s functions.
In the above simulations, I used the
simplGridDBspatial database, because usingsimpleDBcaused an error. However, I’m not sure what the problem was. Will this affect the accuracy of the Green’s function calculation?
If you have points on a grid, it is always better to use SimpleGridDB over a SimpleDB because the interpolation is much faster.
In my simulations, I noticed that multiple files are generated, such as
step07a-leftlowercrust.h5_t2997972000000002.vtk, but they only contain information for theleftlowercrustmaterial. I’m not sure why the other materials are not included.
It looks like you are getting VTK files instead of HDF5 files. You may want to turn off output for the material when computing Green’s functions to reduce disk space and speed up the simulations. Set the output data fields to an empty array [].
Thank you very much for your help. I will give it a try, and if I run into any issues, I’ll reach out to you again. Wishing you all the best!