I am trying to estimate the stress and displacement caused by changes of hydrological loading of a lake. So I am using Pylith to first calculate the surface displacement, following step 18,19 in the manual 2.2.1. The mesh is 1000 km 1000km 100km and the lake is about 100*100. I use the Neumann bc on +z and constrain +x/-x and +y/-y on [0] and [1], and free -z. Also, I specify the normal stress for the region of lake in a SimpleDB. However, the result is unrealistic, where the z-displacement only occurs under the lake but nearby regions have really low z-displacement (please see the attached figure). This very rapid vertical displacement decay does not agree with other modeled results or GPS observations. I am wondering if I missed some steps or any good methods to estimate the surface load like this.
I’m not entirely sure how you are constraining your vertical displacements. Is your bottom boundary free or constrained? If it is free, how are you constraining your vertical displacements? Also, to me it appears that your mesh is quite thin
in the vertical direction, compared to the dimensions of the lake.
Cheers,
Charles
Charles Williams I Geodynamic Modeler
GNS Science I Te
Pῡ Ao
1 Fairway Drive, Avalon 5010, PO Box 30368, Lower Hutt 5040, New Zealand
You are right. I just tried the same load on the different thickness, and I found the modeled z-displacement present different magnitudes. You can see them from attached figures, which illustrate the z-displacements for 100 km, 200 km, and 400 km. However, I am not sure the proper thickness I should choose for this case. For the bottom, I simply constrain the vertical displacement with a fixed z [2]. Is there any suggestion for the bottom bc?
Keep increasing the depth until the solution converges. As a rule of thumb, the boundaries of your domain (sides and bottom) should all be about the same distance from the sources of the deformation. Note that the solution will be quite smooth near the boundaries, so you can use larger cells there. If the displacement and stress field look smooth, then you are likely resolving the solution relatively well. Again, verify convergence of the solution using meshes with different resolutions.
I’m not sure what the proper thickness should be. Is the solution against which you are comparing an elastic half-space solution? I think you will need to consider what you want to measure, and whether the dimensions (and load) of the lake are
large enough that you may need to consider the viscoelastic response to loading (either from the upper mantle or lower crust).
Cheers,
Charles
p.s. I think that fixed-z displacements on the bottom boundary should be OK. One option might be to keep increasing the thickness of your mesh until there is no longer a significant difference in the computed vertical displacements.
Charles Williams I Geodynamic Modeler
GNS Science I Te
Pῡ Ao
1 Fairway Drive, Avalon 5010, PO Box 30368, Lower Hutt 5040, New Zealand
Thank you for the note. I added it to 550 km and the solution converges well and result is at the same magnitude as elastic half-space observation. Also, the material below is set to be vicoelastic.
