Hi Rene,
Thank you very much for your detailed and helpful reply.
First, to answer your questions: for both models, I implement a 2-cm-per-year plate velocity at the top, which is why I added the ‘box with lithosphere’ model in the first case because I can’t prescribe the plate velocity at the top boundary if it is free surface.
For the first case, thanks for pointing out why topography is forced to be zero at x=0 and x=3000km. The reason why I didn’t make lithosphere left and lithosphere right to ‘Additional tangential mesh velocity boundary’ is I implemented the prescribed velocity on those two boundaries.
subsection Boundary velocity model
set Prescribed velocity boundary indicators = bottom: function, lithosphere left: function, lithosphere right : function, front y: function, back y: function, lithosphere front y: function, lithosphere back y: function
subsection Function
set Variable names = x,y,z
set Function constants = vx_top=0.02, vx_bot=0.0, vy_top=0, vy_bot=0, vz_bot=0.0, vz_top=0.0, x0=1000e3,y0=0e3, r_p=200e3, rho_m=3400, alpha=3e-5, g=9.8, delta_T=250, eta0=1e21, gamma=3.1
set Function expression = if(z==0e3, vx_bot, vx_top);0;if(z==0e3,if(((x-x0)^2+(y-y0)^2)<r_p^2, rho_m*alpha*g*r_p^2*delta_T/4/gamma/eta0*exp(-gamma*(((x-x0)^2+(y-y0)^2)/r_p^2))*365*24*60*60,0),vz_top)
end
end
But you’re right. I should also add the lithosphere left and lithosphere right to ‘Additional tangential mesh velocity boundary’.
For the second case, I prescribed velocity on the top boundary and left the left and right boundary ‘open’ (initial lithostatic pressure) so that material can flow in from the left and out of the right boundary. The boundary condition is shown below:
subsection Boundary velocity model
set Prescribed velocity boundary indicators = bottom: function, top: function, front y: function, back y: function
subsection Function
set Variable names = x,y,z
set Function constants = vx_top=0.02, vx_bot=0.0, vy_top=0, vy_bot=0, vz_bot=0.0, vz_top=0.0, x0=1000e3,y0=0e3, r_p=200e3, rho_m=3400, alpha=3e-5, g=9.8, delta_T=250, eta0=1e21, gamma=3.1
set Function expression = if(z==0e3, vx_bot, vx_top);0;if(z==0e3,if(((x-x0)^2+(y-y0)^2)<r_p^2, rho_m*alpha*g*r_p^2*delta_T/4/gamma/eta0*exp(-gamma*(((x-x0)^2+(y-y0)^2)/r_p^2))*365*24*60*60,0),vz_top)
end
end
Yes, it is an ageing plate. The same ageing plate also applies to the first case.
subsection Initial temperature model
set List of model names = adiabatic,function
subsection Adiabatic
# This is the age of the top boundary layer in years.
set Age top boundary layer = 80e6
end
subsection Function
set Variable names = x,y,z
set Function constants = DeltaT=250, r_p=200e3,x0=1000e3,y0=0e3
set Function expression = if(z>=0 && z<=340e3,if(((x-x0)^2+(y-y0)^2)<r_p^2, DeltaT*(1-((x-x0)^2+(y-y0)^2)/r_p^2),0),0)
end
end
subsection Heating model
set List of model names = adiabatic heating
end
This is how the dynamic topography output looks like at 297 Myrs:
The pressure doesn’t change throughout the whole time (300 Myrs) and looks like this:
Thank you very much!
Cheers,
Ziqi
