Boundary traction on spherical geometries with free surface

[mfmweerdesteijn]

When applying boundary traction to the top free surface of a 2D viscoelastic chunk or spherical shell, numerical instabilities occur. Performing the exact same case on a 2D box works very well.

Like in the free surface traction benchmark, a load of length 100 km is applied, increasing linearly from 0 to 100 years. Top boundary is free surface, left and bottom boundaries are free slip, right boundary is "open" using an initial lithostatic pressure. For the chunk the east boundary is free slip and the west boundary is "open" (increasing longitude goes counter clockwise). I show free surface deformation results for the box and chunk model with dimensions 500 x 500 km (box) and 4.5 degree x 500 km (chunk, 4.5 degree is 500 km for an Earth-like radius of 6371 km) and for a much wider case, 20000 x 500 km (box) and 180 degree x 500 km (chunk). Every line in the plot is a 2.5 year time step (from brown to blue).

The chunk free surface deformations seem to be okay-ish from 0 to 4.5 degrees longitude (is 500 km arc distance) compared to the box model, but for the wide case (up to 180 degrees), the free surface heavily oscillates. Also attached are animations of the strain rate through time. For the box models the strain rate increases at the load boundary and they look logical. For the chunk models we see a chaotic pattern which seems to be due to numerical instability.

https://user-images.githubusercontent.com/57097877/134688343-3b23d4a2-83fc-4c0a-9706-6d4f3ce56f7a.mp4 https://user-images.githubusercontent.com/57097877/134688345-e36e6eac-c758-4f3d-822c-a4e86411e4dc.mp4 https://user-images.githubusercontent.com/57097877/134688347-d299e7b5-891c-4293-9901-8fa87534a3fd.mp4 https://user-images.githubusercontent.com/57097877/134688348-9a989ccb-03b4-4f1c-bf12-47c7f72b2326.mp4

Thought: something wrong with coordinate systems in boundary traction application. Test idea: for an increasing chunk radius the chunk curvature decreases more approximates a 'flat' Earth or box model. Result: both for smaller and larger chunk radius the free surface oscillates.

I've also tried with left and right (west and east) closed boundaries, with different mesh refinements, vertical vs normal surface velocity projection for the free surface, but these are not the problem Does anyone have an idea what's going on?:) I'm quite stuck. Any help very much appreciated!

Here are the prm files: case_chunk.txt case_chunk_wide.txt case_box.txt case_box_wide.txt

[tjhei]

I would start with the spherical shell as I am less sure that everything is correct for the chunk (we had many issues on the past). Can you try to create a setup for a shell and share it please? Are you using the free surface stabilization? Does it help to increase the parameter?

[mfmweerdesteijn]

Here is the spherical shell case. The only difference is that I did nullspace removal from the velocity with the angular momentum option, otherwise there's is no convergence. But as you can see the free surface deformation is very oscillating and the strain rates chaotic like in the chunk. Changing the free surface stabilization parameter did not help.

case_shell_nr.txt https://user-images.githubusercontent.com/57097877/134881126-72ea68c0-9bc6-4699-ba61-e927f3694326.mp4

[naliboff]

@mfmweerdesteijn - are there versions of these benchmarks that only consider purely viscous cases? This would be helpful in identifying if the surface of the error is related to the viscoelasticity implementation, surface tractions, or some other feature.

[mfmweerdesteijn]

Okay ja, it's the same prm files but then with a very high elastic shear modulus. I see the same issues in the strain rates and free surface deformation as for the viscoelastic case.

https://user-images.githubusercontent.com/57097877/135057665-2f103a52-e0fd-4935-a05e-16908b6710ee.mp4 https://user-images.githubusercontent.com/57097877/135057670-91efeea6-9413-4b1c-b09b-dbcbf8107ed4.mp4 https://user-images.githubusercontent.com/57097877/135057671-ecc8d958-cd7c-4a65-b9b4-d5e9e0d52252.mp4 https://user-images.githubusercontent.com/57097877/135057672-6fd20d42-ffa2-49aa-b50f-1046aa6c7db6.mp4 https://user-images.githubusercontent.com/57097877/135057674-9f74f34b-261f-4d9f-8fc3-b9241dee6b5e.mp4

[naliboff]

@mfmweerdesteijn - we had a chance to discuss this issue at the user meeting today, and @gassmoeller suggested that may be related to the use of linear mapping when mesh deformation is used. Models with no mesh deformation use a fourth order mapping.

In previous models using a linear mapping with spherical geometries there were pressure oscillations near the surface, which certainly could lead to the type of fluctuations seen here. The small bit of good news is that increasing the resolution reduces the magnitude of these oscillations.

@mfmweerdesteijn, can you try systematically increase the resolution of the current models using the shell geometry, and see if this has any effect?

Another test that may give some insight is to enable mesh deformation, but rather than use the free surface on the upper boundary, have it remain fixed through time. I believe this will still produce the switch from fourth order to linear mapping. We could then see if there are pressure oscillations or other variations in the solution arising from the same time-dependent surface tractions when using a fixed top boundary with or without mesh deformation turned on.

[gassmoeller]

@mfmweerdesteijn, can you try systematically increase the resolution of the current models using the shell geometry, and see if this has any effect?

If that is too expensive, coarsening the resolution should increase the spikes.

Another test that may give some insight is to enable mesh deformation, but rather than use the free surface on the upper boundary, have it remain fixed through time.

Good thought that will separate the effect of the mapping from the effect of the mesh deformation.

[mfmweerdesteijn]

Thank you @naliboff and @gassmoeller that it was brought up in the user meeting and the suggestions here. With the last suggestion, do you mean that the mesh is only allowed to deform tangential to the outer boundary?

[mfmweerdesteijn]

@naliboff @gassmoeller

"Another test that may give some insight is to enable mesh deformation, but rather than use the free surface on the upper boundary, have it remain fixed through time. I believe this will still produce the switch from fourth order to linear mapping. We could then see if there are pressure oscillations or other variations in the solution arising from the same time-dependent surface tractions when using a fixed top boundary with or without mesh deformation turned on."

If I understood it correctly: I used the Mesh deformation subsection, but then set the mesh deformation for the outer boundary to 0 in both directions using 'Boundary function'. I double checked and indeed no surface topography change as intended. In the animation you still see the wild strain rate variations. Does that mean that the problem lies in the mapping and not the deforming mesh?

https://user-images.githubusercontent.com/57097877/135302950-0b279983-3e6b-4bad-af1b-9ffc66bf5fbd.mp4 case_shell_nr_fs.txt

For the increasing or decreasing mesh resolution tests: unfortunately I can't get the solver converging if I do not use adaptive mesh refinements, so I can't systematically increase or decrease resolution to my liking.

[tjhei]

Do we know if the pressure scaling is sensible for these models, @gassmoeller ? That can explain behavior like this.

@mfmweerdesteijn do you run a recent master that includes https://github.com/geodynamics/aspect/pull/3681 ?

[mfmweerdesteijn]

@tjhei Yes I ran it on a version of about a week ago

[naliboff]

@mfmweerdesteijn - Thanks for doing those tests! A few additional questions and associated tests to run:

1. Do the same strain rate fluctuations appear if you apply the same surface traction, but do not have mesh deformation turned on (i.e., the top boundary should still not deform)?
2. Likewise, can you perform similar tests for a cartesian geometry to see if the strain-rate fluctuations appear?
3. Can you also make movies/images of the pressure field?

Hopefully this helps sort out of the underlying issue is the mapping, or something else.

Thanks! John

[gassmoeller]

Do we know if the pressure scaling is sensible for these models, @gassmoeller ? That can explain behavior like this.

I do not immediately see why the pressure scaling should be wrong, the viscosities should be computed correctly, and the representative length scale should not be too different from a realistic model. It's possible because of the thin shell and benchmark, but I would check the other things first.

If I understood it correctly: I used the Mesh deformation subsection, but then set the mesh deformation for the outer boundary to 0 in both directions using 'Boundary function'. I double checked and indeed no surface topography change as intended.

Yes, use the mesh deformation section, but only use the 'boundary function' plugin, not the free surface plugin. Sounds like you did exactly that.

In the animation you still see the wild strain rate variations. Does that mean that the problem lies in the mapping and not the deforming mesh?

It at least means even if the mesh is not deforming you see this problem, I think it is either the elasticity or the mapping.

Okay ja, it's the same prm files but then with a very high elastic shear modulus. I see the same issues in the strain rates and free surface deformation as for the viscoelastic case.

Is it possible to switch of elasticity entirely by setting 'Enable elasticity' to false? It does not need to be a benchmark case, just to see if the oscillations disappear. This would help distinguish between elasticity and mapping.

[mfmweerdesteijn]

1. Do the same strain rate fluctuations appear if you apply the same surface traction, but do not have mesh deformation turned on (i.e., the top boundary should still not deform)?

When I do do that, the solution doesn't converge, already in the first time step.

2. Likewise, can you perform similar tests for a cartesian geometry to see if the strain-rate fluctuations appear?

Turning mesh deformation off in the box model, the solution does converge. Strain rate is similar to with mesh deformation. And you see a pressure increase underneath the surface loading: prm file: case_box_fs2.txt strain rate movie: https://user-images.githubusercontent.com/57097877/135623703-d4cfa767-ec57-430d-9c88-6993a5802938.mp4 nonadiabatic pressure movie: https://user-images.githubusercontent.com/57097877/135623745-f64a87f2-607d-4225-8bc7-b6ab86247973.mp4

3. Can you also make movies/images of the pressure field? nonadiabatic pressure for box, viscoelastic or viscous, and 500x500 or 20000x500 km is exact same as the movie above for the box when mesh deformation is turned off.

For the shell: https://user-images.githubusercontent.com/57097877/135625235-cb8cf568-2abe-4cc5-8d24-a335a2c9ddf3.mp4 (same nonadiabatic pressure pattern for the different setups discussed in this thread)

[mfmweerdesteijn]

Is it possible to switch of elasticity entirely by setting 'Enable elasticity' to false? It does not need to be a benchmark case, just to see if the oscillations disappear. This would help distinguish between elasticity and mapping.

Ja I then used the material model 'simpler model' because for the viscoelastic material model it is required than 'enable elasticity' is set to True. prm file: case_shell_nr_e.txt strain rate movie: https://user-images.githubusercontent.com/57097877/135626042-66004a08-d8ae-4ffd-831c-0d903ee857e3.mp4 nonadiabatic pressure movie: https://user-images.githubusercontent.com/57097877/135626122-1738373c-a410-4ba1-af7f-ceed9ad130c2.mp4 surface deformation plots: The oscillations are nearly identical to when elasticity is enabled for the viscoelastic material model (see plot previously in this thread)

[naliboff]

@mfmweerdesteijn - Thanks you for all of these tests and plots, incredibly helpful!

To make sure I have looked at all pertinent results, is this an accurate summary of the tests to date:

1. Spherical models with surface loading and mesh deformation enabled experience anomalous (unphysical) oscillations, regardless of whether the free surface is enabled or the top boundary does not move.
2. The oscillations are clearly visible in the strain rate , as well as the surface deformation patterns if the top boundary is a free surface.
3. The oscillations are present either with or without elasticity enabled.
4. This oscillations do not appear in cartesian geometry models.

Is all of the above accurate? If so, I think the issue is likely the mapping or another yet to be determined issue.

I'm curious as to why the spherical model with the surface loading, but without mesh deformation is not converging. Was the top boundary set as a zero traction boundary? Would you mind posting the log.txt file for this failed model?

Thanks again for all the work on this! John

[mfmweerdesteijn]

@naliboff All of the above statements are correct!

"I'm curious as to why the spherical model with the surface loading, but without mesh deformation is not converging. Was the top boundary set as a zero traction boundary? Would you mind posting the log.txt file for this failed model?"

I apply surface loading as boundary traction on the top boundary. How can I set that same boundary then as a zero traction boundary? Then there is no surface loading present. Here are the prm file and log.txt: case_shell_nr_fs2.txt log.txt case_shell_nr_fs2_err.txt

What does converge is what I posted before, that the Mesh deformation is enabled but set to zero, but still showing the anomalous oscillations: "If I understood it correctly: I used the Mesh deformation subsection, but then set the mesh deformation for the outer boundary to 0 in both directions using 'Boundary function'. I double checked and indeed no surface topography change as intended. In the animation you still see the wild strain rate variations. Does that mean that the problem lies in the mapping and not the deforming mesh?"

https://user-images.githubusercontent.com/57097877/135832164-1d697db1-95bb-4a54-8d87-6cd25ea9eae4.mp4 case_shell_nr_fs.txt

I'm happy there's a discussion and some momentum going:)

[mfmweerdesteijn]

@gassmoeller @naliboff @tjhei Thanks for the suggestion last meeting to change the mapping in interface.cc of mesh_deformation! I run into a problem that I don't understand, as I think it is quite fundamental FE related, and I was hoping you could have a look at it.

This is the branch: https://github.com/mfmweerdesteijn/aspect/tree/mapping_mesh_deformation I changed the elements order in line 98, and in line 103 I changed to a higher order mapping (from MappingQ1Eulerian to MappingQEulerian).

The error I get is as follows: mapping error.txt When I run: case_shell_nr_map1.txt

[mfmweerdesteijn]

Maybe it has to do with 'spacedim' which is now unspecified and thus spacedim=dim?

From https://www.dealii.org/current/doxygen/deal.II/classTriangulation.html: template<int dim, int spacedim = dim> class Triangulation< dim, spacedim > A triangulation is a collection of cells that, jointly, cover the domain on which one typically wants to solve a partial differential equation. This domain, and the mesh that covers it, represents a dim -dimensional manifold and lives in spacedim spatial dimensions, where dim and spacedim are the template arguments of this class. (If spacedim is not specified, it takes the default value spacedim=dim.)

Thus, for example, an object of type Triangulation<1,1> (or simply Triangulation<1> since spacedim==dim by default) is used to represent and handle the usual one-dimensional triangulation used in the finite element method (so, segments on a straight line). On the other hand, objects such as Triangulation<1,2> or Triangulation<2,3> (that are associated with curves in 2D or surfaces in 3D) are the ones one wants to use in the boundary element method.

Or from https://www.dealii.org/current/doxygen/deal.II/step_34.html#Collocationboundaryelementmethod: Implementation The implementation is rather straight forward. The main point that hasn't been used in any of the previous tutorial programs is that most classes in deal.II are not only templated on the dimension, but in fact on the dimension of the manifold on which we pose the differential equation as well as the dimension of the space into which this manifold is embedded. By default, the second template argument equals the first, meaning for example that we want to solve on a two-dimensional region of two-dimensional space. The triangulation class to use in this case would be Triangulation<2>, which is an equivalent way of writing Triangulation<2,2>.

However, this doesn't have to be so: in the current example, we will for example want to solve on the surface of a sphere, which is a two-dimensional manifold embedded in a three-dimensional space. Consequently, the right class will be Triangulation<2,3>, and correspondingly we will use DoFHandler<2,3> as the DoF handler class and FE_Q<2,3> for finite elements.

[tjhei]

This has nothing to do with spacedim. The default is correct. I am happy to take a look but I won't get to it until the end of the week, probably.

[mfmweerdesteijn]

After a meeting at AGU with @gassmoeller, he found that there is a mistake in source/boundary_traction.function.cc: (line 38-51): if one chooses to apply boundary traction in a spherical coordinate system, then the traction output it also in spherical coordinates. In source/boundary_velocity/function.cc, if a spherical coordinate system is chosen, the coordinates are converted to cartesian (line 53-54). This part is missing in source/boundary_traction.function.cc. The consequence is that an incorrect boundary traction is added to the Stokes system (source/simulator/assemblers/stokes.cc line 932-934).

Here are the radial deformation results, repeated from earlier in this thread, for when we use a spherical coordinate system in the boundary traction subsection (prm file: case_shell_nr_sph.txt )

Here are the radial deformation results for when we use a cartesian coordinate system in the boundary traction subsection (prm file: case_shell_nr_car.txt )

Here is a sketch to better understand the boundary traction application (line 69 in the prm files):

For the cartesian case, the oscillations are still very much present, but of slightly smaller amplitude than the spherical case. The cartesian case shows similar large moving strain rate variations as the paraview animations shown earlier in this thread.

Another thing that was discussed to investigate is not using a custom mesh, but using the default spherical shell, and thicker than what I have now (500 km) because the hyper shell mesh definition works best for certain ratios between shell thickness and circumference. But when I do that, the Stokes solver doesn't reach convergence, already in the first step (2 prm file examples with this problem: case_shell_nr_2.txt and case_shell_nr_3.txt

[anne-glerum]

Hi @mfmweerdesteijn, thanks for the update! Too bad the cartesian way of specifying the boundary traction did not remove the oscillations.

About the full spherical shell models, nothing jumped out to me except maybe the resolution. In the nr_2 case, there is one level of global refinement and one level of adaptive refinement based on the strain rate that is updated every time step. I think it would be better for the test to use a fixed mesh (over time -> set Time steps between mesh refinement = 0) that is more refined. E.g. more global refinement and then some adaptive refinement towards the top. Did you check how the mesh looks in paraview? So these models fail at time 0 already? Is the traction boundary condition still valid for the full sphere?

Does it help when the transition from zero traction to non-zero traction is smooth instead of step-like?

[mfmweerdesteijn]

Thank you for the feedback @anne-glerum!

"Does it help when the transition from zero traction to non-zero traction is smooth instead of step-like?" Is that an option? I now let the boundary traction linearly increase over small time steps, but it is step-like. Is there a smoother/interpolated option?

Indeed with a finer mesh and the cartesian way of specifying the boundary traction I do get convergence, and solutions. I'm still trying out a lot, the mesh refinement seems to play a large role in the solution. Here are a few tests with strain rate animations and surface deformation plots. They are all with constant, but differently refined meshes.

a. mesh resolution of about 100 km case_shell_a.txt b. mesh resolution of about 50 km case_shell_b.txt : no convergence after 62.5 yr c. mesh resolution of about 50 km between +1500 km and -1500 km in phi direction (the load is applied between+100 and -100 km), from 1500 km depth to surface and 100 km elsewhere case_shell_c.txt d. mesh resolution of about 50 km between +1500 km and -1500 km in phi direction (the load is applied between+100 and -100 km), from CMB to surface and 100 km elsewhere case_shell_d.txt : no convergence after 70.0 yr e. mesh resolution of about 12.5 km between +300 km and -300 km in phi direction, from 300 km depth to surface, mesh resolution of about 50 km between +1500 km and -1500 km in phi direction, from 1500 km depth to surface (apart from more refined region) and 100 km elsewhere

strain rate animations (note different scales): a. https://user-images.githubusercontent.com/57097877/151994969-7516f6c7-bf9d-440f-8d4c-85e0c29191a1.mp4 b. https://user-images.githubusercontent.com/57097877/151995019-50e9b715-c6ee-4834-9c10-a7b1b36780e5.mp4 c. https://user-images.githubusercontent.com/57097877/151995047-797d670d-cd15-48f3-82d0-11fead708184.mp4 d. https://user-images.githubusercontent.com/57097877/151995072-2c7d0e2b-a071-4687-ae30-a366980c0f65.mp4 e. https://user-images.githubusercontent.com/57097877/151995098-23566a77-5345-4b41-9da5-39be38fcccd7.mp4

surface deformation plots: a.

b.

c.

d.

e.

Seen in the animations: when the mesh is constant throughout the domain, we see strain rates at the surface where the boundary traction is applied. When the mesh has different refinements throughout the domain, we see very large strain rates (much larger than the ones at the surface from the boundary traction) at the corners of these subdomains. Seen in the zoomed in deformation plots: The more refined the mesh, the less of a spikey, and more rounded deformation is visible, which is good. But with these more refined domains we still see unexpected deformation patterns, likely to do with these large strain rates at the corners of the more refined regions. I do not understand these strain rates. This is what the deformation should look like (modeled in a box geometry, and not necessarily the right values, just an idea for the shape):

[anne-glerum]

Hi @mfmweerdesteijn, good to see the progress!

What is the resolution in km for the box geometry? For a load that is applied from -100 to + 100 km, I do think you need a resolution more on the order of model e or smaller. But that's probably hard for a uniform mesh. Maybe you can still use the normal spherical shell, but with a smaller opening angle than 360?

The strain rate artefacts at the mesh refinement corners I've seen before in chunk models.

Concerning the step-like traction boundary conditions, through the function expression one can make a smooth transition from zero to max traction, using for example a htan() function.

[mfmweerdesteijn]

Hello @anne-glerum!

in the box it's 3.25 km near the ice load (12.5 km for spherical case e here). I think I can use a spherical geometry with smaller opening angle, but to have such a refined mesh throughout the whole domain, and then also in 3D (considering the GIA application), it becomes computationally very expensive, especially compared to other codes, so that's not very appealing.

Okay so you've come across a similar issue before in chunk geometries. A good thing to discuss next week!

Thank you for the function tip. I am wondering, doesn't it then still only evaluate at the numerical time step and be step like? Im not sure I understand this approach.

[mfmweerdesteijn]

Hi @naliboff, could these strain rates at the corners of the higher mesh resolution have something to do with elasticity in spherical geometries? I have never seen this in box geometries with different levels of refinement. Also, I assume these strain rate anomalies are not encountered when modeling mantle convection and plumes etc, what aspect is build for. What do you think?

[mfmweerdesteijn]

@anne-glerum @naliboff

I redid case c, but with only viscous material (elasticity not enabled). case c viscoelastic (VE): case_shell_c.txt (repeated from above) case c viscous (V): case_shell_c_v.txt

Strain rate animations (note, they have the same scale): VE: https://user-images.githubusercontent.com/57097877/152960718-89a4e602-b6d2-49af-aac2-259e9bd49971.mp4 V: https://user-images.githubusercontent.com/57097877/152960793-4fe74844-9dbf-4091-b0a3-6a925b082cd8.mp4

The strain rate anomalies at the boundaries of the refined mesh are present in both cases, but much smaller in the V case than the VE case. (You can see the strain rate at the surface from the boundary traction with this scale as well for both cases)

Vertical surface deformation through time: VE:

V:

The strain rate anomalies for the Viscous case are large enough still and affect the surface deformation...

[mfmweerdesteijn]

Update:

1. Do we get good deformation results (no jumps) when mesh is radially refined (such that we do not get the strain rate anomalies which occur for refined regions with lateral boundaries)? case b: global mesh resolution of about 50 km case_shell_b.txt case f: mesh resolution of 50 km up to 300 km depth, 100 km mesh resolution elsewhere (radially refined mesh) case_shell_f.txt

surface deformation plots through time: b: f:

strain rate animations: b: https://user-images.githubusercontent.com/57097877/154036755-52e46d77-aecd-487c-a192-44044e5ba8b4.mp4 f: https://user-images.githubusercontent.com/57097877/154036783-85ee5111-5e5d-453a-8b29-7e67ec21582d.mp4

Surface deformation result is identical, although strain rates are different: for the radially refined mesh (f) we see a large strain rate at the interface between the different refined regions, but apparently this does not affect the surface deformation solution, as it did with the strain rate anomalies at the lateral boundaries for more refined regions.

2. Do the strain rate anomalies at the refined regions lateral boundaries get less significant when we increase the surface loading? Or, do the strain rate anomalies scale with the surface loading? case c: mesh resolution of 50 km between +1500 km and -1500 km in phi direction (load is applied between+100 and -100 km), up to 1500 km depth, 100 km mesh resolution elsewhere case_shell_c.txt case g: same as case c but with 100x increased surface load case_shell_g.txt

surface deformation plots through time: c: g:

strain rate animations: c: https://user-images.githubusercontent.com/57097877/154039615-4220a669-19a0-4430-a214-2f7563445a9c.mp4 g: https://user-images.githubusercontent.com/57097877/154039640-2e06b234-2bc8-4fb2-afc3-233b6360b750.mp4

In case c the applied ice height is 1 km and in case g 100 km. In the surface deformation plots we see that the solution for case g is more smooth (much smaller jumps due to the strain rate anomalies). For case c the max. deformation is 14.5 m and for case g 335 m. We use linear viscoelasticity so we would expect a 100x difference in surface deformation as we surface load is x100 difference. However, 335 m/100 = 3.35 m. Thus, the strain rate anomalies have less of an impact for larger surface loading; the strain rate anomalies do not linearly change with the surface loading.

3. Does the solution with a refined region with lateral boundaries improve when we use a higher order mapping for the spherical geometry? case c: mesh resolution of 50 km between +1500 km and -1500 km in phi direction (load is applied between+100 and -100 km), up to 1500 km depth, 100 km mesh resolution elsewhere case_shell_c.txt case h: same as c but performed on PR #4379 in which spherical geometry is approximated with higher order mapping case_shell_h.txt

surface deformation plots through time: c: h:

strain rate animations: c: https://user-images.githubusercontent.com/57097877/154040941-d3644dcc-e226-4104-824a-f4639d743b40.mp4 h: https://user-images.githubusercontent.com/57097877/154041455-3f845d8b-4b75-4614-99b2-dc8df4433acb.mp4

Running the same prm file on PR #4379 with the higher order mapping leads to odd sharp jumps in the surface deformation solution.

[mfmweerdesteijn]

I started comparing the surface deformation solution for a 2D box and shell geometry. Because of the high resolution layer required in the spherical geometry (to avoid the strain rate anomalies at boundaries of change in mesh resolution), I went with a chunk, to save computational resources.

Surface load is a 100 km length line (of ice) in 2D, which linearly increases from 0 to 1 km ice height over 100 yr. The lateral boundary at the load is a free-slip boundary (axisymmetric case). We use a wide geometry, to minimize the effect of the far boundary. Mesh resolution between box and chunk is very similar. Upper 200 km is ~3.2 km resolution and below is ~50 km resolution. Width, or arc distance, of the models is 3000 km, depth is 500 km.

case 1 box and chunk: far boundary is also free-slip. case 2 box and chunk: far boundary is 'open' by applying an initial lithostatic pressure.

CASE 1 box prm file: case_box_1w.txt chunk prm file: case_chunk_1w.txt

box surface deformation:

chunk surface deformation:

We see that the surface deformation magnitude for the chunk is slightly smaller than for the box.

CASE 2 box prm file: case_box_1wo.txt chunk prm file: case_chunk_1wo.txt

box surface deformation:

chunk surface deformation:

Again, we see that the surface deformation magnitude for the chunk is slightly smaller than for the box. But there is also uplift in the middle field, and more in the far field next to the 'open' boundary, which shouldn't be there.

chunk, far lateral boundary, strain rate animation: https://user-images.githubusercontent.com/57097877/154697692-f8b51f58-8c63-4ee4-944c-8bad18d40187.mp4

An unexplained strain rate anomaly which appears when using a far 'open' boundary. This anomaly does not appear in the box geometry. Bug in initial lithostatic pressure implementation for spherical models?

[mfmweerdesteijn]

@anne-glerum @gassmoeller @naliboff Previously I showed that the surface deformation for the 2D chunk gets close to that of the 2D box with free-slip lateral boundaries, and a radially or vertically refined mesh (no lateral mesh region boundaries). Why is the surface deformation not equal? This could be dependent on the chunk radius.

Shown here: surface deformation for box model, and for chunk model with increasing radius from left to right: The maximum surface deformation (at t=100 yr, end of sim., dark blue line) for the chunk of R=0.5xRe is 67.1% of the max surface deformation of that of the box. for R=Re this value is 90.8% and for R=2xRe it is 97.4%. With a larger radius the chunk better approximates a "flat earth" and thus box geometry. This is great. But, for R=4xRe the surface deformation is all over the place, and we see large strain rate anomalies at the bottom boundary: https://user-images.githubusercontent.com/57097877/155147483-38c498d7-956b-411b-8e3f-4473fea7b03e.mp4 I don't understand why the solution breaks down for a larger radius.

Another observation from before: the surface deformation for the 2D chunk with an 'open' boundary (initial lithostatic pressure) at the far side of the load shows a large strain rate anomaly at that open boundary at the bottom, resulting in an uplift feature at the boundary. How does this uplift feature behave for different chunk radii? We see that this uplift feature gets smaller with an increase in radius. At R=2xRe there is still some light oscillation occurring, and again at R=4xRe the solution breaks down with large strain rate anomalies at the bottom boundary of the chunk. But as we see this uplift feature getting smaller with an increasing radius (thus more flat model), could there be a bug in the implementation of the initial lithostatic pressure in spherical geometries?