Weird velocity vectors

ASPECT
1

Hello everyone,

I am trying to create a compressional environment. The deformation is initiated through consistent inward horizontal velocities at the model sides, without imposing any velocity limitations in the vertical directions. Outflow at the model base equalizes lateral inflow. The inward velocity is 1 cm yr-1. The bottom boundary has a free slip, the top boundary is a free surface. The balance calculations look okay to me. However, my velocity vectors look weird, they rapidly change direction in each timestep. I tried to increase free surface stabilization theta to 1, enlarged the model box to triple sizes, decreased CFL number to 0.1 and max time step to 1e4, and decreased velocity down to 0.05 cm yr-1. They did not affect anything.
outputs.pdf (1.1 MB)

I would appreciate any suggestions to fix this issue. Thanks in advance.

Best,
Acelya

type or paste cosubsection Initial temperature model
  set Model name = function
  subsection Function
    set Variable names = x,y
    set Function constants = h=660e3, ts1=273, ts2=441, ts3=723, ts4=1623, \
                             A1=1.8e-6, A2=0.25e-6, \
                             k1=2.5, k2=2.5, k3=2.5, \
                             qs1=0.0415, qs2=0.031, qs3=0.0125
    set Function expression = if( (h-y)<=15.e3, \
                                  ts1 + (qs1/k1)*(h-y) - (A1*(h-y)*(h-y))/(2.0*k1), \
                                 if( (h-y)>15.e3 && (h-y)<=40.e3, \
                                      ts2 + (qs2/k2)*(h-y-15.e3) - (A2*(h-y-15.e3)*(h-y-15.e3))/(2.0*k2), \
                                       if( (h-y)>40.e3 && (h-y)<=220.e3,\
                                         ts3 + (qs3/k3)*(h-y-40.e3), ts4)))
  end
end

type or paste code subsection Boundary velocity model
  set Tangential velocity boundary indicators =  bottom
  set Prescribed velocity boundary indicators = left x: function, right x:function

  subsection Function
    set Function constants  = cm=100.0 #*3600.0*24.0*365.25
    set Function expression = if(x<1650e3 && z<370.e3, -(230/430)/cm, \
                                  if(x<1650e3 && z>360.e3, 1.0/cm, \
                                    if( x>1650e3 && z<370.e3, (230/430)/cm, \
                                        if(x>1650e3 && z>360.e3, -1.0/cm, 0) \
                                      ))); \
                                

    Constant viscosity prefactors                         = 0.1,                                0.1                                    
    set Prefactors for dislocation creep          = background:6.52e-15,  C_1:6.52e-15
    set Stress exponents for dislocation creep    = background:3.5  ,     C_1:3.5
    set Activation energies for dislocation creep = background:530.e3  ,  C_1:530.e3 
    set Activation volumes for dislocation creep  = background:18.e-6,    C_1:18.e-6
   
    set Prefactors for diffusion creep            = background:6.52e-14  , C_1:6.52e-14                     
    set Grain size exponents for diffusion creep  = background:0.        , C_1:0    
    set Activation energies for diffusion creep   = background:530.e3      C_1:530.e3  
    set Activation volumes for diffusion creep    = background:3.e-6   ,   C_1:3.e-6    
    set Cohesions                                 = background:0.e6    ,   C_1:0.e6             
    set Angles of internal friction               = background:15.0      , C_1:15.0           
    set Grain size = 5.e-3
    

  end
end
2

Dear Acelya,

Can you please attach your prm file? It is easier to assess what is wrong when we can see what you’ve done.

Best wishes,
Bob

3

Hi Bob,

I attached .prm file too.

original.prm (10.0 KB)

4

Thanks for the prm, that is very helpful!

I’ll defer to someone with more experience of the viscoplastic model for detailed advice (@anne-glerum or @jbnaliboff may be able to help, though they are both very busy right now). The following are a few less-expert observations.

It is clear that your simulation is not conserving mass - this is particularly obvious in the second figure in your outputs.pdf.

Thoughts:

  • You might want to simplify your model to pin down where the problem is. You have a viscoplastic model with a free surface - it would be good to know whether the problem arises from one of those two choices, or the combination of the two.
  • Do the results change if you reduce the solver tolerance? ASPECT’s default values are significantly smaller than what you have chosen in the prm.
  • You have set the rock cohesion to 0. If you set it to a more normal rock-like value, does that have an effect? Zero cohesion means that your rock at the surface is globally yielding, and I don’t know whether we have tested this.

Best wishes,
Bob

5

Dear Acelya,

I agree with Bob, try to simplify your model setup, in particular the boundary condition at the upper surface. From what you prescribe at the side boundaries, one would expect more outflow at the top, and it does seem strange that you get these large changes between time steps.

It looks like you also have hillslope diffusion at the upper surface, so you could try to see if the instability is still there if you switch that off, or if you simply leave the boundary open rather than using a free surface, just to figure out which of the features in your model causes the problem.

Juliane

6

Dear @bobmyhill and @jdannberg ,

Thank you for your time and suggestions. I will try them!

When I use free slip on the surface, inflowing material looks fine, but outflowing still shows instabilities (maybe this is normal?).

I tried to simplified my input file and removed the crust (usually, crust has cohesion value of 1e6). In both cases, I have the same problem.

Best,
Acelya

7

Dear Acelya,

When I use free slip on the surface, inflowing material looks fine, but outflowing still shows instabilities
(maybe this is normal?).

“Flickering” in deformation style between timesteps is normally an indication that either the solver is struggling or the simulation is poorly resolving the features of interest. It isn’t normal.

One other thing you should check carefully, especially if you change the top boundary to be free slip, is your boundary velocity problem:

set Function expression = if(x<1650e3 && z<370.e3, -(230/430)/cm,
if(x<1650e3 && z>360.e3, 1.0/cm,
if( x>1650e3 && z<370.e3, (230/430)/cm,
if(x>1650e3 && z>360.e3, -1.0/cm, 0)
))); \

  • Clearly define behaviour at every depth. You can use <= to indicate “less than or equal to” to avoid using different bounds (instead of using 360 and 370 km, you should choose a single depth and apply it)
  • Sharp steps in velocity are a bad idea. Much better to have a continuous velocity gradient (e.g., adding a region in between your inflow and outflowing regions where the inward velocity varies linearly with depth).
  • If you intend to have no net outflow (e.g. by imposing free slip on the top surface), you need to ensure that your constraints can satisfy mass conservation. At the moment, you have a net flux of -230/430/100*370000 = -1979 m^2 / year through the bottom left of your domain, and
    1/100*(660000 - 370000) = 2900 m^2 / year through the top left of your domain. This net influx is inconsistent with mass conservation in a closed domain. This is what @jdannberg was using to make this statement:

From what you prescribe at the side boundaries, one would expect more outflow at the top

Best wishes,
Bob

8

Dear Bob,

I tried to create a 10 km gap between inflowing and outflowing materials. But I just realised that I made the wrong calculation!
Based on your suggestions, I reconducted the model and the crazy velocity vectors are gone!
Thank you again Bob.

Best,
Acelya