Dear everyone,
I am new to PyLith (4.2.1) and have been testing some simple examples related to the problems I hope to study, starting with box-3d.
The issue I am running into is the following. If I set up a compression model with shortening only in the x direction, while the other two directions are fixed, the model thickens during compression. However, because Uz = 0 at the bottom boundary, this thickening seems to be expressed only as surface uplift, which may be overestimated because there is no root growth. My understanding is that this may not be very consistent with an Airy isostatic compensation.
So my question is: does PyLith have a boundary-condition setup that is more suitable for this kind of problem, where compressional thickening can be accommodated by both upward uplift and downward thickening? Whether any other boundary-condition options in PyLith might provide a more realistic basal constraint, since the support at the base would presumably evolve as the crust thickens, rather than remaining constant. This makes me think that prescribing a single fixed force at the bottom may also be problematic.
PyLith does not have special boundary conditions for the type of simulation you describe. It has Neumann and Dirichlet boundary conditions. For elasticity, these correspond to traction and fixed displacement boundary conditions.
To model the type of behavior you describe, there are multiple options. You could use a traction boundary condition on the bottom boundary, but it would need to be in equilibrium with the rest of the domain. Another approach is to make the model deep enough so that a Dirichlet boundary condition on the bottom makes sense; this is usually accompanied by variations in the material properties.
Thank you for the clear explanation.
For Option 1, as the model deforms over time in a viscoelastic simulation, the bottom traction needed for isostatic balance would need to change with the displacement — essentially requiring a displacement-dependent boundary condition rather than a fixed or time-dependent traction. I think this makes it difficult to implement in PyLith.
For Option 2: My understanding is that by extending the model deep enough (e.g., 400–600 km) and assigning low-viscosity material (η ~ 10^17–10^18 Pa·s) near the bottom, the weak layer would effectively decouple the Dirichlet boundary from the near-surface deformation and mimic isostatic compensation.
One thing I am curious about: in such a setup under horizontal compression, would the model produce downward velocities in the mid-depth or near-bottom region (i.e., the thickening expressed as crustal root deepening, rather than purely as surface uplift)? This would be a key indicator that the weak layer is indeed providing isostatic-like behavior.
Option 1: You are correct. Some codes implement “Winkler forces” (effectively springs on the bottom boundary) to approximate this behavior. We have not implemented this in PyLith, because we prefer Option 2.
Option 2: The details of the behavior depend on the details of the variation of the material properties with depth. The desired behavior can be region dependent; in some regions, you might expect more thickening (uplift at the surface) than in others. I think the best approach is to capture the main features you would expect in a self-gravitating 3D ellipsoid (Earth).
Thank you, that makes sense. One thing worth noting is that Earth’s viscosity structure still carries significant uncertainties, and these uncertainties may affect the modeled uplift patterns differently depending on the choice of boundary conditions. Anyway, I plan to run some tests to explore it.
