We are currently conducting a simulation study on the surface deformation of a viscoelastic medium under the combined influence of gravitational and tectonic loading. During our investigation, we referenced the “subduction (3D)” module (PyLith 2.2.1) in relevant case studies and executed the following procedures: First, we ran step08a to obtain the surface deformation induced by gravity (Figure 1). Subsequently, we used generate_initial_stress.py to introduce gravitational forces into the material model and performed step08c. However, upon analyzing the simulation outcomes, we observed perplexing phenomena: at time zero, the domain’s surface deformation was negligible (Figure 2), but by the second time step, the surface subsidence abruptly became substantial (Figure 3), and it continued to escalate with time.
We question the accuracy of these results. In a simulation framework that incorporates both viscoelasticity and gravity, how should the simulation be correctly configured? Additionally, when applying gravitational forces to the material, we have encountered uncertainties: is it sufficient to load only the six stress components, or is it necessary to simultaneously load the six stress components, total-strain components, and viscous-strain components?
Figure 1
Figure 2
Figure 3
If you want to transfer the state from one simulation to another, then you need to include all state variables. For a viscoelastic material, this includes the stress tensor, viscous strain, and total strain.
There are several possible approaches for combining viscoelasticity and gravity. It is often challenging to construct realistic simulations because the Earth’s stress state is the result of millions of years of deformation. We generally know the geometry at the current time but not the state of the system. You can try to reconstruct the deformation to achieve an estimate of the current state of the system, or you can generate a rough estimate based on simple assumptions.
A completely relaxed state as an initial state for further modeling
You can run a simulation with gravity for a sufficient number of time steps so that the shear stresses are negligible. A difficulty with this approach is that the deformation resulting from the relaxation may substantial so that the geometry at the end of the simulation does not match the current geometry that you want to use for further modeling.
Elastic stress state for further modeling
You can run a simulation with gravity and an elastic material to estimate the stress state as a starting point for viscoelastic deformation. This is often a reasonable approach when the stress state is close to the lithostatic overburden pressure (no shear stress). In this case, you only transfer the stress field from the simulation with gravity and an elastic material to the simulation with viscoelasticity. If the stress state has significant shear stress, then you should expect relaxation of the shear stress in the simulation with viscoelasticity.
Stress state from incompressible elasticity for further modeling
Using PyLith v3 and later (we recommend the current release), you can use the incompressible elasticity formulation to obtain a stress state wth minimal deformation. This approach works well if you want to estimate the stress state associated with topography and the shallow materials are elastic. If you need features available in PyLith v2 for simulation of further deformation, then you can use the current release to compute the initial state and generate the input needed for a simulation with PyLith v2; we recognize this is klunky and we are working hard to reimplement all of the features in v2 in the current version.
Viscoelasticity and finite strain
It is challenging to obtain a stable numerical solution for problems with viscoelasticity, gravity, and finite strain. A stable solution requires high accuracy, especially for the stress field. This is particularly challenging in PyLith v2, which is limited to a basis order of 1; this is discussed in PyLith v2 in examples/2d/gravity (this suite of examples illustrates several approaches). PyLith v3 and later do not yet have finite strain, so this approach is not available.
