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
A you have discovered, producing desired convection patterns with boundary conditions is not always straightforward
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?
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.
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
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
below, I upload pictures of the 1st, 2nd, 3rd, 4th and 16th timeframe of the resulting velocity field. As you can see, the plume gets dragged to the right, contrary to the initial conditions specified.
thanks for your suggestions. In the code snippet I now introduced a velocity field that varies linearly between -0.045 and 0 between y=2800km and y=2600km by specifying
set Function expression = if(y>2.8e6, -0.045, 0) + if(y<=2.8e6 & y>2.6e6, -(y - 2.6e6)*0.045/0.2e6, 0) ; 0
However, this didn’t result in the plume getting dragged to the left, as we want it. Afterwards, I shifted the linear decay of the leftward translation a bit by changing the numbers 2.8e6 and 2.6e6 above to 2.5e6 and 2.0e6 respectively. Still, this didn’t improve the results. So all in all, unfortunately, I still cannot explain why the plume in the frame keeps getting dragged to the right (see the pictures I deployed above) although I specify left-oriented translations as initial conditions.
As a reference point, it would be really helpful to see the velocity field produced by the BCs above if there is no plume (i.e., don’t prescribe a thermal anomaly). Could you run a model that does not contain any thermal perturbations to confirm that the prescribed velocity field is as expected?
Also, the velocity values you are prescribing (4.5 cm/year) may be low relative to the velocity field produced by thermal anomaly. Can you show a plot of the velocity field magnitude when you only use tangential velocity BCs, and the flow field is driven purely by variations in buoyancy?
This should help provide some constraints on what prescribed velocity values will be needed to deflect the plume.
and sorry for the late reply.
I looked at your input file, and there one thing I would like to point out:
You don’t prescribe a top boundary condition. Is that intended? It means that your mantle is open at the surface and material can flow in and out. If you do what John suggested (having a lithosphere with a prescribed velocity), then you should also specify the velocities at the top boundary, or at least make to top impermeable.
Why does the plume in your model gets dragged to the right? I think that makes sense: You prescribe flow at the surface of the models that goes towards the left, but the side and bottom boundaries are otherwise closed. So to conserve mass, the bottom of the mantle has to move to the right (because the top is moving to the left). Once your plume approaches the surface, it might be dragged towards the left.
If you want to see the influence of mantle flow on the plume, I would prescribe a velocity gradient over the whole side boundary (for example, being zero at the bottom of the box and your maximum value at the top of the box). This should lead to your plume being dragged in that direction (as long as the top boundary is closed).
In addition, I wanted to link to two papers where we have used boundary conditions from Bernhard Steinberger’s global models to deflect plumes in regional convection models. They only look at the upper mantle and not the whole mantle though. But these papers both have data packages as supplement that contain the code and input files to reproduce the models (they are a few years old, so might not work with the current ASPECT version, but could be useful examples).
I think closing the top boundary by setting a tangential velocity boundary indicator for the top actually solved the problem of the plume getting dragged into the wrong direction. Currently, I’m satisfied by the results of the simple model and will get back to you if there come up any questions in the following.
after some trying around, I think there’s one question left open for me concerning the results that I obtained so far. As I pointed out in the last post, by closing the top boundary, the desired effect on the plume is now visible. However, I’m getting artifacts after about 100 timesteps in the calculations that lead to abortion of the program shortly after. Below, I attached my current code and one picture, showing the artifacts. Do you have any idea why these come up and how to prevent this?
Thank you very much in advance,
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 = 5e6
set Output directory = NGP-vm5c
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
#vertical velocity gradient, one plume dragged to the left
subsection Boundary velocity model
set Tangential velocity boundary indicators = top
set Prescribed velocity boundary indicators = left x: function , right x: function
subsection Function
set Variable names = x,y
set Function expression = -0.35*y/3e6 ; 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 Boundary heat flux model
# set Fixed heat flux boundary indicators = 0,3
#subsection Function
#set Function expression = 0.2
#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
I think there are at least two different problems in the model:
You have an inflow velocity at the right boundary of the model, but you do not prescribe a temperature (or composition) at that boundary. That leads to the huge temperatures you see at that boundary.
Your resolution might not be high enough. I would try to rerun the model with a higher refinement level and test if that improves the temperature oscillations you see around the plume.
