initial lithostatic pressure in spherical geometries

[mfmweerdesteijn]

When applying initial lithostatic pressure profiles on far boundaries (to make them 'open') in spherical geometries (here 2D chunk), large strain rate anomalies occur at these boundaries, resulting in unwanted free surface deformation (if used in combination with free surface, as is the case here). Here are plots of the free surface through time (colors) for a linearly increasing load with time at the center of the domain. Ignore the plot titles.

1: boundary traction and closed lateral boundaries

2: boundary traction and open lateral boundaries

2: strain rates for case 2

3: no boundary traction, only open lateral boundaries

3: the latter but zoomed in

The initial lithostatic pressure causes strain rate anomalies at the applied boundaries, in spherical geometries. These strain rate anomalies can be of significant size w.r.t. the boundary traction. Is this a mapping problem? Any ideas?

[naliboff]

Thanks @mfmweerdesteijn! So, it looks like a small amount of deformation is occurring at the boundaries even when no surface traction is applied (i.e., the model should be static). Is this interpretation correct?

If so, can you add an image of the strain rates for the "open boundaries, no surface tractions" case? If there are strain anomalies, can you add an image of the pressure field and adjust the color scale to highlight small differences at a specific depth?

Thanks!

[mfmweerdesteijn]

Indeed, there are deformations even when no surface traction is applied, so coming from the initial lithostatic pressure at the lateral boundaries applied.

open boundaries, no surface tractions strain rate images: START (t=0)

END (t=100 yr)

I played around with the pressure field but didn't manage to visualize small differences. But I'm doing a test now with equal resolution throughout the domain, because it seems like some of the high strain rates are located near the border of the mesh resolution.

[mfmweerdesteijn]

No boundary traction, only initial lithostatic pressure on both lateral boundaries, equal mesh resolution throughout domain: again strain rate anomaly, at bottom boundary, resulting in unwanted free surface deformation.

[anne-glerum]

Hi @mfmweerdesteijn , I went through the initial_lithostatic_pressure plugin, and I didn't see anything obvious why the computed pressure wouldn't be correct in the spherical case. The chunk geometry does have its issues though. Do you see the same anomalies in a quarter spherical shell? Also, there are some checks in the initial_lithostatic_pressure plugin that are only checked when run in Debug mode. Have you done that?

[mfmweerdesteijn]

@anne-glerum thanks for looking into it! I think I run in Debug mode yes, let me check.

When I use a spherical shell with 90 degrees opening angle and with initial lithostatic pressure on the lateral boundaries, I get the error: "After displacement of the mesh, this function can no longer be used to determine whether a point lies in the domain or not." This is from the point_is_in_domain function. But I don't see why it would work for other geometries (box and chunk) with mesh deformation enabled. I then changed the conditions of the AssertThrow in the spherical shell to the same conditions as in the chunk.

from:

template <int dim>
bool
SphericalShell<dim>::point_is_in_domain(const Point<dim> &point) const
{
  AssertThrow(!this->get_parameters().mesh_deformation_enabled ||
              this->get_timestep_number() == 0,
              ExcMessage("After displacement of the mesh, this function can no longer be used to determine whether a point lies in the domain or not."));

to (based on conditions in the point_is_in_domain function from the chunk geometry):

template <int dim>
bool
SphericalShell<dim>::point_is_in_domain(const Point<dim> &point) const
{
  AssertThrow(!this->get_parameters().mesh_deformation_enabled ||
              this->get_timestep_number() == 0 ||
              this->get_timestep_number() == numbers::invalid_unsigned_int,
              ExcMessage("After displacement of the mesh, this function can no longer be used to determine whether a point lies in the domain or not."));

These are btw the conditions in the box geometry:

template <int dim>
bool
Box<dim>::point_is_in_domain(const Point<dim> &point) const
{
  AssertThrow(!this->get_parameters().mesh_deformation_enabled ||
              this->simulator_is_past_initialization() == false,
              ExcMessage("After displacement of the free surface, this function can no longer be used to determine whether a point lies in the domain or not."));

Why isn't the code consistent?

With this change, the model runs, and unfortunately I still get strain rate anomalies at the bottom boundary at the lateral boundaries, and at a few random along the bottom boundary:

surface deformation (ignore plot title):

strain rate (one capture of the 90 degrees shell with strain rate anomaly in the corner at the lateral boundary, and an anomaly randomly placed along the bottom boundary):

[anne-glerum]

Hi @mfmweerdesteijn thank you for the additional info.

Why isn't the code consistent?

These functions have ben improved on over time and have apparently not been updated everywhere. The last implementation (for the box) is the most correct. It would be great if you can put in a PR or make an issue to remind us to unify the functions.

The strain rate anomalies in the bottom corners: Couldn't they derive from the difference in velocity directions? The bottom boundary is tangential, but velocities along the open boundaries can have tangential as well as normal components (I think). So at the corner, there could be a change in direction, leading to a higher strain rate. Can you check what the velocity vectors look like there?

Also, is there flow into or out of the domain (on both open boundaries)? Then the question is, why is the computed pressure not good enough and induces flow? What are the velocity magnitudes?

The surface deviations are very small (less than 1 m right), especially compared to the model resolution (what is the resolution?). Maybe they become smaller with higher resolution. By default a 1000 integration points are used to compute the lithostatic profile. How many elements are there in the depth direction?

You can also analytically compute the lithostatic pressure at the bottom of your model and compare it with what is computed in the plugin (with a print statement). When the density is constant, this is just densitygravitydepth.

an anomaly randomly placed along the bottom boundary

Are there more along the bottom boundary? What happens to the velocity there?