Strange pattern on viscosity field

Hi everyone,

I am working on a continental collision model where an oceanic plate subducts beneath a continent while dragging another continental piece which eventually results in a continental collision. The model seems to be running fine but I am getting this grid-like pattern in the viscosity field after the initial timestep which makes me believe something is wrong. See the figure for clarification.

Example viscosity field:

I am using Geodynamic World Builder to create geometry. 6 km thick oceanic crust on the subducting plate is modified to act as a lubricant layer (I decreased the viscosity of this layer artificially to create sufficient lubrication). You can find .prm and .wb files below. Other than the examples below; I have tried to switch off the free surface, I have tried models with just diffusion creep and composite flow laws but the results are similar to the example screenshots.

.prm file:

#### World Builder parameters which can be used by ASPECT ######## World Builder parameters which can be used by ASPECT ####
set World builder file = sub2col_base_2PlatDip45.wb

set Dimension = 2
set Start time                             = 0
set End time                               = 50e6
set Pressure normalization                 = no
set Maximum time step                      = 1e4
set Maximum first time step                = 100
set Maximum relative increase in time step = 20
set Use years in output instead of seconds = true

set Nonlinear solver scheme                = iterated Advection and Stokes
set Max nonlinear iterations               = 30
set Nonlinear solver tolerance             = 1e-3
set CFL number                             = 0.1
subsection Solver parameters
  subsection Stokes solver parameters
    set Linear solver tolerance                = 1e-6
    set Number of cheap Stokes solver steps    = 100

set Output directory = sub2col_base_v31

subsection Initial temperature model
    set Model name = world builder

subsection Boundary temperature model
  set Fixed temperature boundary indicators = bottom, top
  set List of model names = box
  subsection Box
    set Bottom temperature = 1600
    set Top temperature    =  273

subsection Initial composition model
    set Model name = world builder

#### parameters needed by ASPECT when using compositions ####
subsection Compositional fields
   set Number of fields = 5
   set Names of fields = Oceanic_Plate, Mantle, Plateau_Crust, Plateau_Mantle, Oceanic_Crust

subsection Geometry model
  set Model name = box
  subsection Box
    set X extent = 2640e3
    set Y extent = 660e3
#    set Z extent = 660e3
    set X repetitions = 8
    set Y repetitions = 2
#    set X periodic = true

subsection Mesh refinement
  set Initial adaptive refinement        = 0
  set Initial global refinement          = 5

##Free Surface
subsection Mesh deformation
  set Mesh deformation boundary indicators = top: free surface
  subsection Free surface
    set Free surface stabilization theta = 0.6
    set Surface velocity projection = vertical

##Boundary classifications
subsection Boundary velocity model
  set Tangential velocity boundary indicators = bottom, left, right

subsection Material model
  set Model name = visco plastic

  subsection Visco Plastic
    set Viscous flow law = dislocation

    # Reference temperature and viscosity
    set Reference temperature = 293
    set Reference viscosity = 1e22

    # The reference strain-rate is used on the first non-linear iteration
    # of the first time step when the velocity has not been determined yet.
    set Reference strain rate = 1.e-16
    set Minimum strain rate = 1.e-20

    # Limit the viscosity with minimum and maximum values
    set Minimum viscosity = 1e18
    set Maximum viscosity = 1e23

    # Material properties
    set Thermal diffusivities = 1.333333e-6
    set Heat capacities       =        750.

                                       #Oceanic_Plate, Mantle, Plateau_Crust, Plateau_Mantle, Oceanic_Crust
    set Densities             =        3300, 3300, 3300,  2900, 3300, 3300
    set Thermal expansivities =        2e-5

    # Harmonic viscosity averaging
    set Viscosity averaging scheme = harmonic

    set Prefactors for dislocation creep          =   6.52e-16,   6.52e-16,  6.52e-16,  8.57e-28,  6.52e-16, 5.5e-18,
    set Stress exponents for dislocation creep    =      3.5,     3.5,       3.5,       4.0,       3.5,      3
    set Activation energies for dislocation creep =   530.e3,     530.e3,    530.e3,    223.e3,    530.e3,   345.e3
    set Activation volumes for dislocation creep  =   18.e-6,     18.e-6,    18.e-6,     0.,       18.e-6,   38.e-6

    # Plasticity parameters
    set Angles of internal friction =   5.7392, 5.7392, 5.7392, 5.7392, 5.7392, 5.7392
    set Cohesions                   =   1e6
    set Maximum yield stress        =   4e8


subsection Gravity model
  set Model name = vertical

# Post processing
subsection Postprocess
  set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization, topography, dynamic topography
  subsection Visualization
    set List of output variables = density, viscosity, strain rate, error indicator
    set Output format = vtu
    set Time between graphical output = 5e5
  subsection Topography
    set Output to file = true
    set Time between text output = 5e5

.wb file:

  "cross section":[[0,330e3],[2640e3,330e3]],
  "coordinate system":{"model":"cartesian"},

    {"model":"oceanic plate", "name":"Subducting Plate", "max depth":100e3,
     "temperature models":
       {"model":"plate model", "max depth":100e3, "spreading velocity":0.01,
       "ridge coordinates":[[0e3,-1],[0,2000e3]]}
     "composition models":[
     {"model":"uniform", "compositions":[4], "min depth":0, "max depth":6e3},
     {"model":"uniform", "compositions":[0], "min depth":6e3,"max depth":100e3}]

    {"model":"continental plate", "name":"Plateau_2", "coordinates":[[800e3,0],[800e3,660e3],[1200e3,660e3],[1200e3,0]],
     "temperature models":[{"model":"linear", "max depth":120e3, "bottom temperature":1600}],
     "composition models":[
     {"model":"uniform", "compositions":[2], "min depth":0, "max depth":20e3},
     {"model":"uniform", "compositions":[3], "min depth":20e3,"max depth":120e3}]

    {"model":"continental plate", "name":"Plateau_1", "coordinates":[[2640e3,0],[2640e3,660e3],[1660e3,660e3],[1660e3,0]],
     "temperature models":[{"model":"linear", "max depth":120e3, "bottom temperature":1600}],
     "composition models":[
     {"model":"uniform", "compositions":[2], "min depth":0, "max depth":20e3},
     {"model":"uniform", "compositions":[3], "min depth":20e3,"max depth":120e3}]

     {"model":"mantle layer", "name":"upper mantle", "min depth":100e3, "max depth":660e3,
        "temperature models":[{"model":"uniform", "temperature":1600}],
        "composition models":[{"model":"uniform", "compositions":[1]}]

     {"model":"subducting plate", "name":"Slab",
      "dip point":[1e7,-1e7],
          {"length":250e3, "thickness":[100e3], "angle":[0,30]},
          {"length":250e3, "thickness":[100e3], "angle":[30,45]}],
      "temperature models":[{"model":"linear", "max distance slab top":100e3, "top temperature":273, "bottom temperature":1600}],
      "composition models":[
      {"model":"uniform", "compositions":[4], "min distance slab top":0, "max distance slab top":6e3},
      {"model":"uniform", "compositions":[0], "min distance slab top":6e3,"max distance slab top":100e3}]    }



That viscosity pattern certainly doesn’t look right. Let me start with asking a few questions/suggestions.

  1. What is your minimum grid size? Is the grid refinement changing a lot each refinement step?
  2. Could you take a look at your temperature and compositional fields? Do they have large under or overshoots?
  3. Could you check whether you actually reach the nonlinear solver tolerance you set of 1e3?
  4. A nonlinear solver tolerance of 1e-3 is not that strict. Could you try setting it to 1e-4 (and make sure you reach that) and see if that solves or changes the issue?

I hope this helps with finding the issue!

Hi Menno,

Thanks for your reply.

  1. The minimum grid size is around 5 km and I only use global refinement.
  2. I can’t see any issue with the temperature and compositional fields.
  3. For nonlinear solver tolerance of 1e-3 model seems to be converging well.
  4. I have done some further testing for nonlinear solver tolerance of 1e-4 and the models do not converge well even with increased maximum number of nonlinear iterations and decreased maximum timesteps. I would appreciate any suggestions to make models converge better.

I also did some further testing with different viscosity averaging schemes and I noticed that using geometric viscosity averaging helps with the problem, but I still do not understand why the harmonic would not work.


Sorry for the late reply.

  1. Having a uniform grid may prevent some potential issues, so it is good that you started with that.
  2. good
  3. good
  4. That is a bit worrying. How many nonlinear iterations did you try? At about what nonlinear residual does it get stuck, or does it keep going down slowly?
  5. Yes, geometric viscosity averaging will probably create less viscosity contrasts in the cold slab. Harmonic averaging favors low viscosities, while geometric averaging favors higher viscosities. This will probably lead to smoother viscosities in those regions. Could you show a figure of that?

Does this problem occur from the start or does it develop later in the model? If it develops later, what are the dynamics of the model while it occurs?

I do notice that you have very strong viscosity contrast in some regions in the model, are these desired?