Cannot achieve desired effect of mantle flow on plume through boundary indicators

Hello there,

I’m currently using version 2.2.0 of ASPECT and I’m trying to model a plume whose path is driven by mantle flow (see Steinberger et al. (2019): Yellowstone Plume Conduit Tilt Caused by Large-Scale Mantle Flow, esp. figure 3 and 4). In order to achieve this I’d like to implement a velocity field split in two, namely a left-oriented velocity field in the lower part of the figure and a right-oriented velocity field in the upper part of the figure, using the function Prescribed Velocity Boundary Indicators. However, in the simulation, the plume immediately gets dragged to the right and the effect of the lower boundary condition just described is not really visible anymore and, in addition, the setup leads to multiple convection cells I didn’t aim for. Below you can find the code I used.

Do you have any ideas what I can do in order to see the effect of both velocity field directions on the plume, both in the left and in the right direction? Thank you in advance for your answers!

Dina

set Dimension = 2

set Start time                             = 0
set End time                               = 2e7
set Use years in output instead of seconds = true
#set Maximum time step			   = 1e7
set Output directory                       = NGP-vm5a

set Use conduction timestep = true

subsection Geometry model
  set Model name = box

  subsection Box
    set X repetitions = 4
    set Y repetitions = 1
    set X extent = 12e6
    set Y extent = 3e6
  end
end

##################### Boundary conditions



subsection Boundary velocity model
	set Prescribed velocity boundary indicators = left x: function , right x: function
	subsection Function
		set Variable names = x,y
		set Function constants = cm=0.1
		set Function expression =if(y<1.5e6, -0.45*cm, 0.45*cm); 0
	end
	set Zero velocity boundary indicators = bottom
end


subsection Boundary temperature model
  set Fixed temperature boundary indicators   = bottom, top
  set Model name = box
  subsection Box
    set Top temperature = 1600
    set Bottom temperature = 4000
  end
end

subsection Initial temperature model
  set Model name = function
  subsection Function
    set Variable names      = x,z
    set Function expression = 3100 - z*0.5e-3 + if (x>9.5e6 & x<10e6 & z<0.5e6,2000,0)
  end
end

subsection Material model
  set Model name = multicomponent

  subsection Multicomponent
    set Densities = 3300,3600
    set Viscosities = 1e21,1e22
    set Viscosity averaging scheme = harmonic
  end
end

## Use compositional field to simulate a lower mantle of different composition
subsection Compositional fields
  set Number of fields = 1
end

subsection Initial composition model
  set Model name =  function
  subsection Function
    set Variable names = x,z
    set Function expression = 1
  end
end

subsection Gravity model
  set Model name = vertical
  subsection Vertical
    set Magnitude = 10.0
  end
end

subsection Heating model
  set List of model names = shear heating
end

subsection Mesh refinement
  set Initial adaptive refinement        = 0 
  set Initial global refinement          = 5
  set Time steps between mesh refinement = 10
  set Strategy = temperature,density
end

subsection Postprocess

  set List of postprocessors = visualization

  subsection Visualization

    set Number of grouped files       = 0
    set Output format                 = vtu
    set List of output variables      = density,viscosity,heating
    set Time between graphical output = 1e5
    set Interpolate output = true
  end
end

Hi Dina,

Welcome and thank you for posting to the forum!

A you have discovered, producing desired convection patterns with boundary conditions is not always straightforward :slight_smile:

A general question: How did you arrive at the current set of boundary conditions? Are they based off the methods used in the Steinberger paper? As there are a number of studies in the literature examining the effects of mantle flow on plumes, I think the most efficient approach is to try and replicate one of the studies in ASPECT and then diagnose why differences are arising.

If this is the approach you are already taking, can you provide the exact BC conditions you are trying to replicate in ASPECT?

Thanks!
John

Hi John,

thank you for your answer. My current boundary conditions are not yet based on those in the Steinberger paper or any other source out there as I just tried to implement the general idea of the paper in a simpler fashion. However, as you said, I think I need more background information on possible boundary conditions in order to obtain a desired pattern, if even possible.

Regards, Dina

Hi Dina,

Sounds good and please let us know as soon we can answer additional questions.

As discussed, in this case we can likely provide the most help with setting up models to reproduce previous results.

Thank you again for posting the question to the forum!

Cheers,
John

Hello John,

I further reduced the complexity of the model so that I’m currently trying to model only one velocity direction as a boundary condition, namely a left-oriented velocity field over the whole area. The goal is now that the plume gets dragged to the left by the velocity field. In the code snippet below you’ll find this minimal example. However, the exact opposite as expected happens and the plume gets dragged to the right. Do you have any ideas why this happens and what can be altered to prevent this effect?

Additionally, do you have any other papers than the Steinberger et al. paper in mind that can be used as an orientation and that contain examples that serve as a good basis for implementation?

# Box model with adiabatics and different composition and all that good stuff

set Dimension = 2

set Start time                             = 0
set End time                               = 2e7
set Use years in output instead of seconds = true
set Maximum time step			   = 5e6
set Output directory                       = NGP-vm

set Use conduction timestep = true

subsection Geometry model
  set Model name = box

  subsection Box
    set X repetitions = 4
    set Y repetitions = 1
    set X extent = 12e6
    set Y extent = 3e6
  end
end

##################### Boundary conditions

subsection Boundary velocity model
	set Tangential velocity boundary indicators = bottom
	set Prescribed velocity boundary indicators = left x: function , right x: function
	subsection Function
		set Variable names = x,y
		set Function constants = cm=0.1
		set Function expression = -0.45*cm ; 0
	end
end

subsection Boundary temperature model
  set Fixed temperature boundary indicators   = bottom, top
  set Model name = box
  subsection Box
    set Top temperature = 1600
    set Bottom temperature = 4000
  end
end

subsection Initial temperature model
  set Model name = function
  subsection Function
    set Variable names      = x,z
    set Function expression = 3100 - z*0.5e-3 + if (x>9.5e6 & x<10e6 & z<0.5e6,2000,0)
  end
end

subsection Material model
  set Model name = multicomponent

  subsection Multicomponent
    set Densities = 3300,3600
    set Viscosities = 1e21,1e22
    set Viscosity averaging scheme = harmonic
  end
end

## Use compositional field to simulate a lower mantle of different composition
subsection Compositional fields
  set Number of fields = 1
end

subsection Initial composition model
  set Model name =  function
  subsection Function
    set Variable names = x,z
    set Function expression = 1
  end
end

subsection Gravity model
  set Model name = vertical
  subsection Vertical
    set Magnitude = 10.0
  end
end

subsection Heating model
  set List of model names = shear heating
end

subsection Mesh refinement
  set Initial adaptive refinement        = 0 
  set Initial global refinement          = 5
  set Time steps between mesh refinement = 10
  set Strategy = temperature,density
end

subsection Postprocess

  set List of postprocessors = visualization

  subsection Visualization

    set Number of grouped files       = 0
    set Output format                 = vtu
    set List of output variables      = density,viscosity,heating
    set Time between graphical output = 1e5
    set Interpolate output = true
  end
end

Dina: could you attach a picture of the velocity field you get?

Separately, this isn’t right:

		set Function constants = cm=0.1

:slight_smile:

Best
W.

Hi Dina,

I think the boundary conditions in your model will simply produce a leftward translation as they extend from the top of the model to the base.

Rather, I think what you want is horizontal velocities along the sides that vary in magnitude as you move from the lithosphere down into the convecting mantle. This will produce shear that deflects the plume.

For example, you could have the leftward (e.g., negative horizontal velocities) in the upper 100-200 km, which then linearly gradate to a horizontal velocity of 0 over the next 100-200 km. This would just require using the model height (y) in the function.

Does this make sense?

Others in the community - can you think of classic papers that have performed this type of simulation? I know there must be some recent examples, but this is not my focus area :slight_smile:

Cheers,
John