Negative surface elevation throughout the whole surface in spherical shell & free surface deformation

ASPECT
1

Hi all,

I am trying to use free surface in 2D spherical shell with ice rheology.

However, I currently see two points that look like numerical errors.

  1. Surface elevation ends up with negative values throughout the top surface.

Here are temperature in domain and surface elevation along the boundaries at 5 Myr.

free_T
free_T4000×1700 110 KB

And below is the velocity magnitude and direction of velocity. Arrows only show direction, scale is not included.

free_v
free_v4000×1700 367 KB

Surface velocity is not large, in order of 1e-5 m/yr, but it shows radially negative values and sinking.

However, I don’t think it is right evolution of surface, as there should be surface uplift balancing sinking.

  1. There are velocity anomalies where mesh size changes by adaptive mesh refinement.

This is the same image from previous velocity image but now it shows zoomed at lower left corner with mesh in red.

free_v2
free_v24000×1700 1.86 MB

There are spots of velocity values with highs and lows, and it happens where mesh size changes.

I think these spots affects convection development, as values are near order of convection velocity.

I have tested same parameters in box geometry and I don’t see these issues in box models.

Also, null space removal methods were tried but it fails to converge at the very beginning.

Could I get some help to remove these two issues? I attached the prm file just in case.

europa_chunk_45_ref_newE_60_iso_vp_lowvis2_free_tests.prm (7.5 KB)

Thanks for your help.

Cheers,

Hyunseong

2

Hi @hyunseong96 ,

Thanks for sharing the detailed explanation and your data. From your description, it sounds like the issues with surface elevation and velocity anomalies near mesh refinement boundaries could be related to mesh deformation handling.

I recently worked on a fix that adds better support for higher-order finite element degrees during mesh deformation, which might help resolve these numerical artifacts. You can find the code in Fix mesh deformation for GMG by adding support for higher-order FE degrees

Please let me know if you encounter any problems or need help applying the patch.

Best,
Ninghui

3

Hi @tiannh7 ,

Thanks for the reply! I will update the ASPECT version and try that!

Meanwhile, I have few questions from your fix.

  1. Does this only work for GMG solver?
  2. Does the code mean that default velocity discretization for free surface deformation was Q1 before and now is Q2?

Thank you again for the help!

Cheers,

Hyunseong

4

Hi Hyunseong,

  1. Does this only work for GMG solver?
    The fix applies to both GMG and AMG solvers. However, there is some additional support specifically for higher-order finite element degrees related to mesh deformation implemented for GMG.
  1. Does the code mean that default velocity discretization for free surface deformation was Q1 before and now is Q2?
    Yes, regarding the velocity discretization for free surface deformation, the previous default was Q1, and now it uses Q2 elements. This helps better capture curved geometries with higher accuracy.
5

@hyunseong96 - Thanks for posting this question to the forum. In addition to the suggestions from @tiannh7, would you mind adding mass flux statistics to the list of post processors?

As you noted, it looks like the entire surface is moving down, but for your model setup (free-slip on all other sides) this would fundamentally violate mass conservation. I would be curious to see what this post processor actually calculates for the mass flux across the boundary.

John

6

Thanks both @tiannh7 @jbnaliboff !

I have updated to the latest version of ASPECT and used GMG solver. However, I kept seeing the minus topography in the end.

However, seeing the statistics about outward mass flux and total mass does tell why I am having leaps in mass flux. It happens when I change AMR to coarser level. First on is total mass graph and the second graph shows mass flux at top and bottom boundaries with time.

iso_60_free_mass_c
iso_60_free_mass_c5000×3000 250 KB

iso_60_free_mass_flux_c
iso_60_free_mass_flux_c5000×3000 400 KB

There is change of mesh refinement level at 0.3, 0.5 and 1 Myr. I changed maximum refinement level from 5 to 2 at each time. I lowered this value so I can have smaller DOFs and larger time step size, which is regulated by CFL number and mesh size.

I guess I should work more on AMR. Meanwhile, if you know any other ways I should try, please let me know! :slight_smile:

Cheers,

Hyunseong

7

@hyunseong96 From your input file, I see that you are using a free surface at the top, and tangential velocity boundary conditions at the left, right, and bottom. So there really shouldn’t be any mass flux through the bottom, but nonetheless you are seeing a mass flux through that boundary. That’s a contradiction that is relatively obvious and that one can debug by itself. My suggestion would be to just keep making the problem simpler and simpler:

  • Does the contradiction remain if you switch to a constant viscosity model?
  • Does the contradiction remain if you switch off mesh refinement?
  • Does the contradiction remain if you <… some other simplification …>

The idea is that you have identified something that shouldn’t happen, but does happen. Now let’s reduce the complexity of the problem so that perhaps we can isolate the underlying issue. Does that make sense?

I will note that it’s conceivable that what you see is just numerical noise. Over the course of your simulation, you lose about 1% of mass, but you’re doing at least 100,000 time steps if I read the input file correctly (and 7000 time steps for the time shown in your figure). So in each time step, you lose around 1e-6 of the mass. That shouldn’t happen, but it’s a pretty small amount. My first simplification for the list above would therefore be to substantially tighten the nonlinear and linear solver tolerances, and see whether that has an effect. For example, change the nonlinear solver tolerance from 1e-4 to 1e-6 and check whether that affects the mass balance. Do the same for the linear solver tolerance.

Best
WB