Dynamic Topography seems to be uneffected by Nullspace removal

ASPECT
1

Hello everyone,

I’m doing 3d simulations where I evaluate the dynamic topography. I got results I was not expecting and thought that the unremoved Nullspace net rotation was the problem. But after I removed the net rotation the resulting dynamic topography is the same. Why is it not effected?Because the flow pattern looks drastically different.

The original problem is that I got a result for dynamic topography for a model where there is no density contrast (or temperature) only viscosity. (I know this is very unrealistic but I wanted to see what this would result in.) As I understand it if there is no lateral density variation the dynamic topography should be zero. The following skript gives a non zero dynamic topography that has quite a complex pattern even though the model is very simple. When a do the same with a density anomaly instead of viscosity the result looks reasonable.

Thanks in advance, Maja

Here is my Skript:

A simple setup for convection in a 3d shell. See the

manual for more information.

set Dimension = 3
set Use years in output instead of seconds = true
set End time = 1000
set Output directory = output-mu_23_const_rm

subsection Material model
set Model name = multicomponent

subsection Multicomponent
set Densities = 4400
set Viscosities = 1e22, 1e23
set Viscosity averaging scheme = harmonic
end
end

subsection Compositional fields
set Number of fields = 1
end

subsection Initial composition model
set Model name = function
subsection Function
set Coordinate system = spherical
set Variable names =r,phi,theta
set Function expression = if(r>3480000 & r< 4480000 & phi > 0 & phi < pi/2 & theta>pi/4 & theta<3pi/4 || r>3480000 & r< 4480000 & phi > pi & phi < 3pi/2 & theta>pi/4 & theta<3*pi/4,1,0)
end
end

subsection Initial temperature model
set Model name = function
subsection Function
set Coordinate system = spherical
set Variable names =r,phi,theta
set Function expression = if(r>3480000 & r< 4480000 & phi > 0 & phi < pi/2 & theta>pi/4 & theta<3pi/4 || r>3480000 & r< 4480000 & phi > pi & phi < 3pi/2 & theta>pi/4 & theta<3*pi/4,1500,1500) #erst llsvp, dann mantel
end
end

subsection Geometry model
set Model name = spherical shell

subsection Spherical shell
set Inner radius = 3480000
set Outer radius = 6371000
set Custom mesh subdivision = list of radial values
set List of radial values = 3586033.3333333335, 3682066.6666666665, 3778100.0, 3874133.3333333335, 3970166.6666666665, 4066200.0, 4162233.333333333, 4258266.666666667, 4354300.0, 4450333.333333333, 4546366.666666666, 4642400.0, 4738433.333333333, 4834466.666666666, 4930500.0, 5026533.333333333, 5122566.666666666, 5218600.0, 5314633.333333333, 5410666.666666666, 5506700.0, 5602733.333333333, 5698766.666666666, 5794800.0, 5890833.333333333, 5986866.666666666, 6082900.0, 6178933.333333333, 6274966.666666666
set Initial lateral refinement = 2
end
end

subsection Boundary velocity model
set Tangential velocity boundary indicators = top, bottom
end

subsection Boundary temperature model
set Fixed temperature boundary indicators = top, bottom
set List of model names = spherical constant

subsection Spherical constant
set Inner temperature = 1500
set Outer temperature = 1500
end
end

subsection Nullspace removal
set Remove nullspace = net rotation
end

subsection Gravity model
set Model name = radial constant
end

subsection Mesh refinement
set Initial global refinement = 1
set Initial adaptive refinement = 1
set Strategy = temperature
set Time steps between mesh refinement = 15
end

subsection Postprocess
set List of postprocessors = visualization, dynamic topography, geoid
subsection Visualization
set Output format = vtu
set Time between graphical output = 1e6
set Number of grouped files = 1
set List of output variables = density,viscosity
end

subsection Depth average
set Time between graphical output = 1.5e6
set Output format = vtu
end
end

subsection Checkpointing
set Steps between checkpoint = 50
end

2

Hi @MajaZimmer: Would you mind showing us the results you get, so that we don’t have to run the simulation ourselves?
Best
W.

3

mu_23

I am not sure what exactly you want me to show you. This is the map of dynamic topography from the simulation without nullspace removal. The one with nullspace removal looks exactly the same but I can upload only one file. I dont know where this pattern comes from and I dont understand why the nullspace removal does not have an effect. I can see in paraview that the nullspace removal worked. Could it be that the removal of the nullspace takes places after the dynamic topography is calculated?
I hope this helps. Do you want me to show you anything else?

Best, Maja

4

Hi Maja,

I believe what you’re seeing are the results of numerical instabilities that form at some of the grid boundaries. There are two reasons that these give you such a large signal in your specific setup. (1) You’re running the model at very low resolution (global refinement of 1). The dynamic topography (and flow) calculation will not be accurate at this resolution. A good way to check that the DT calculation is more accurate is to look at the global mean (which should be close to 0). (2) You’re not prescribing any lateral density perturbation, which I believe enhances the instabilities. Once you prescribe density perturbations these will kick off flow and reduce this grid related pattern. These density perturbations can be very long wavelength (e.g. degree 2) if you want to keep the model simple.

Best,
Jacky

5

Hi Maja,
I think Jacky already described the origin of the anomalies, it is likely that they are a mesh artifacts and show you exactly in which places the mesh conforms more or less to an ideal sphere (they should change/decrease if you change the resolution).

I just wanted to mention one additional point about the nullspace removal: While nullspace removal is important for models that have a nullspace (i.e. models that have no side boundaries and tangential velocities at top and bottom), I do not think it would ever affect the dynamic topography you see. That is because the nullspace removal removes the part of the solution that lies in the nullspace of the Stokes operator, which means the part of the solution vector that will not affect the solution of the equation. In practice this is mostly a solid-body rotation (in a spherical geometry) or solid-body translation (in a box geometry). So in other words the part of the solution that gets changed close to the surface by nullspace removal is only flow that is tangential to the surface. But only flow that is not tangential to the surface will cause dynamic topography. I am simplifying a bit here, lateral pressure gradients may contribute as well, but they are not changed by the nullspace removal either. Thus it wouldnt surprise me if nullspace removal would never change the dynamic topography you see in the output of instantaneous models. In models that run over time it would of course affect where material is transported to in later timesteps, so in those models it could matter.

Best,
Rene

6

Hi Maja,

The pure rotation from nullspace removal is essentially a degree-1 toroidal flow in 3D spherical geometry. Therefore, it’s a horizontal flow along the surface, but the radial velocity gradient (radial stress) is used to calculate dynamic topography. What Rene described clearly detailed this point.

In addition, I’d like to add one more point about your rheology. From your prm file, you divide the viscosity domain to two areas with the values 1e22 and 1e23, respectively (multicomponent material model). If you compare the boundaries between these two compositional areas defined in your prm file and your plot of dynamic topography. You can see that the bizarre large dynamic topography lies along these boundaries. This actually does not surprise me because the abrupt horizontal jump of the viscosity will indeed induce unreasonably large dynamic topography wiggles across the boundary. This is another numerical phenomenon in calculation of the dynamic topography in finite elements across the lateral viscosity jump. I also see this in my models. You can try to increase your resolution to see whether these dynamic topography wiggles can be reduced (or narrowed down).

Can you describe more on your purpose to turn off temperature anomaly while turning on the lateral viscosity division of your domain? People commonly start from a simple model with a simple temperature anomaly pattern (like a single wavelength) and a constant (or radial) viscosity field.

Best,
Shangxin

7

Hi Jacky,

thank you very much for your answer. The tip for the mean of the dynamic topography is greatly appreciated. I increased the global refinement to 2 and the mean DT is now closer to zero and the overall amplitude of the DT is decreased to around 10 m. When I include density variations the outcome looks much more reasonable especially with increased resolution.
Thanks again!

Best,
Maja

8

Hi Rene,

Thank you for the explanation of the Nullspace removal. It makes a lot of sense now.
I increased the resolution and the amplitude of the anomalies got smaller as you predicted.
Thanks again!

Best,
Maja

9

Hi Shangxin,

Thank you for your answer. I increased the resolution and the dynamic topography decreased.
There is not really a purpose for this model :smiley: I was just curious what the results would look like and then I was confused by the pattern. But it makes much sense now that it is only a numerical problem. In my other models I use density and temperature anomalies with different viscosity structures for the mantle like you descibed.
Thanks again for your reply!

Best,
Maja