I would like to investigate dynamic topography at the CMB in spherical geometry. This is typically done with a free-slip boundary, where topography is calculated based on stresses. Yet, it would be interesting to use a “free base”, similar to the models which study surface topography.
When I had a look at the manual (section Free surface), it sounds like you could apply that to any boundary of your domain. However, all set-ups that I have tried ran into a solver convergence error sooner or later, no matter what choice of “CFL number” and “free surface stabilization theta” I used.
I then tried the cookbook “free-surface.prm” and switched the free surface to the bottom, but again I end up with the convergence error after 3 time steps with a huge residual (~1e20).
Has anyone ever tested using the free surface at any other boundary than the top of the domain? And is this convergence error related to the fact that you cannot specify any material properties (e.g. density, viscosity) for the fluid outside your domain? I guess ASPECT intrinsically uses the values for air, which might cause problems for the CMB.
Broadly, a ‘free surface’ implies that the boundary is zero stress and in ASPECT the pressure at the free surface boundary is zero (e.g. no overlying material).
In the case of assuming a free surface at the CMB, unrealistic velocities and lack of convergence is likely due to material moving rapidly towards the model center. In other words, without a density contrast or imposed pressure at the CMB, I would predict that material would just move radially towards the model center and the magnitude of this motion would swamp out anything else.
My initial thought for how to model development of CMB topography is to use a “sticky molten-iron approach”. This is similar to the sticky-air approach, where a low-viscosity (ex: 10^18 Pa s verse 10^26 Pa s) and low-density (air or water) layer is placed above the lithosphere. The boundary between these layers thus behaves as a “quasi free-surface”.
This method is used extensively in the geodynamics community and nearly exclusively by codes using finite-difference methods. For more details, see the following paper:
In the case of the CMB, you would impose a low-viscosity but high-density layer beneath the mantle.
In fact, it looks Taras Gerya may have tried a similar approach for modelling propogation of iron diapirs within silicate materials: https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/ggge.20078
Hope this helps and please keep us updated with any additional questions or new results you would like to share!
Looking forward to hearing how that approach works out.
A few notes of caution regarding the “stick-air” approach:
Having a very large viscosity contrast (ex: 10^8) will likely require quite a bit more work by the solvers, even for a linear rheology. I suggest starting with a smaller viscosity contrast (10^2 or 10^4) and then systematically increasing the magnitude from there.
In order for the boundary between the mantle and “sticky iron” to behave in a realistic manner, you need to have sufficient resolution across the boundary and be far enough away from the basal boundary.
AMR will likely be a big help here. FYI, this issues are all discussed in the first paper listed from my first response.
Hi Björn,
John’s suggestion (sticky-iron) is one possible solution to the problem, however at least theoretically you are right, the moving free surface approach should also work for the CMB. I do not think anyone has done this particular problem, and it might be some work to get it to run, it is hard to determine if it would be more or less work than the sticky-iron approach.
The trick to using the free surface approach for the bottom boundary is to exactly revisit what the free surface approach is doing:
Prescribe a stress boundary condition at one boundary: For the surface you simply remove all boundary indicators from velocity or traction boundary conditions, which as John wrote implies the stress is zero. For the CMB this is not the right approach. You would need to prescribe a realistic pressure to get realistic velocities through the boundary, and you can use the prescribed traction boundary conditions for that. The initial lithostatic pressure plugin might just work out of the box, but you might need to play around with the function plugin to prescribe a correct pressure. This is the first step that you can test without any of the steps below (more stable).
If the Stokes solution from the above step is realistic (no oscillations, realistic flow) these velocities can be used to actually deform the boundary. Add the bottom boundary to the free surface indicator, which will cause the boundary to be deformed by the computed velocities.
The machinery in ASPECT should handle everything else correctly if 1. and 2. worked (i.e. adjusting the rest of the mesh to avoid distorted cells, correct velocities to still solve the correct advection equation).
Depending on how much time you want to spend on this, you might give at least 1. a try.
Setting up this sticky metal layer is indeed a bit tricky, as you already pointed out. A large viscosity contrast is very likely to cause a solver convergence error, even for rather high mesh resolutions. The same holds true when trying a very high thermal conductivity in order to make that layer more or less isothermal.
To get some first test running, I had to play with the AMR strategy, because you need a good resolution throughout that sticky metal layer to avoid solver instabilities. I talked to Fabio Crameri about all that since he has some experience with the sticky air approach.
thank you for suggesting a strategy how it could be possible to make the free surface approach applicable to the bottom.
I will have a closer look at the respective plugins and probably test your step 1 to check whether it is worth pursuing this further. It would be really nice to have a “free base” approach working properly, but I have to see how much time I can actually spend on this.
I will keep you updated on any advance in this direction.
I tried your suggestion and step 1 works reasonably well with the initial lithostatic pressure plugin.
However, once you try to use the free surface, you receive an error message from the simulator/helper_function telling you that you have assigned multiple velocity or traction boundary conditions to one boundary. So it seems that the current implementation of the free surface in ASPECT automatically sets the pressure (or more precisely the traction) to 0 and does not allow to change that.
Is there an easy way around that, e.g. telling the helper function to ignore that, or would this cause problems somewhere else in the code?
Hi Björn,
I think this is just a consequence of nobody ever tried this. I think this error message is incorrect, it is usually not the right thing to combine a free surface with velocity or traction boundary condition (if the free surface were at the top), but for your case it is the right thing to do. So I would just comment out the variable parameters.free_surface_boundary_indicators from the assignment to boundary_indicator_lists in the first lines of the function Simulator<dim>::check_consistency_of_boundary_conditions in helper_functions.cc. That should silence the error.
I switched off the error message and the code runs. However, it obviously ignores the Free surface boundary indicators since the bottom boundary remains an open boundary and the elements are not deformed. So it seems to me that, when trying to combine the traction boundary and the free surface, the traction boundary has higher priority.
Do you have any suggestion how to solve that problem? For example, is there any way to modify simulator/free_surface.cc in a way that you can define a pressure at the boundary?
You are right, I had a look through the free surface code, to see where traction boundaries are removed. You will need to comment out the following code block in source/simulator/free_surface.cc:FreeSurfaceHandler<dim>::make_constraints():
// Zero out the displacement for the traction boundaries
// if the boundary is not in the set of tangential mesh boundaries
for (std::map<types::boundary_id, std::pair<std::string, std::string> >::const_iterator p = sim.parameters.prescribed_traction_boundary_indicators.begin();
p != sim.parameters.prescribed_traction_boundary_indicators.end(); ++p)
{
if (tangential_mesh_boundary_indicators.find(p->first) == tangential_mesh_boundary_indicators.end())
{
VectorTools::interpolate_boundary_values (free_surface_dof_handler, p->first,
ZeroFunction<dim>(dim), mesh_displacement_constraints);
}
}
That should allow traction boundaries to be moved around with the velocity. Let me know how that goes.
I have been experimenting with the free base, but there is still a problem I am not able to solve.
The changes you mentioned allow to have both traction boundary conditions and a free surface. However, setting the correct traction at the boundary is challenging. The initial lithostatic pressure plugin does not do the job since it cannot account for the pressure variations induces by the topography. Therefore, I describe the traction as a function (based on the theory of the sticky air approach)
P_0 + drho ⋅ g ⋅ (3481e3 - R)
with drho being the density difference between mantle and core and R being the spherical coordinate of the function plugin. P_0 is an offset pressure that I introduced while trying to stabilize the algorithm.
Since I use the Boussinesq approximation, P_0 = 0 is the obvious choice. However, this set-up causes the whole CMB to be uplifted over time, which causes numerical instabilities to occur after a certain while and eventually leads to a breakdown of the solver. Varying P_0 from positive to negative values does not solve the issue (positive values cause even faster breakdown and completely unrealistic topography oscillations).
I though that the missing pressure normalization at the surface might be a problem since I could get negative pressures within the domain, so I switched off the respective check in core.cc. However, this does not have any effect, except for the pressure field being positive everywhere.
Do you have any idea how to solve that problem? Does it relate back to how the pressure is calculated and used for the Boussinesq approximation within ASPECT?
The changes you mentioned allow to have both traction boundary conditions and a free surface.
Ok, that is great to know.
However, setting the correct traction at the boundary is challenging. …
Yes that is the tricky part. Your approach sounds right though. It is hard to tell from here what might go wrong, but have you checked if the prescribed traction at the boundary leads to the same pressure at the bottom boundary as in a model without the free bottom boundary?
Since I use the Boussinesq approximation, P_0 = 0 is the obvious choice.
This I do not quite understand. Should not P_0 be set to the pressure that is usually expected at the CMB, and the variable part would then model the pressure variations due to the deformation of the boundary? Note that even when ASPECT uses the Boussinesq approximation we solve for the full pressure (unless you significantly modified it to only use the dynamic pressure). I would run a model without open bottom boundary (e.g. tangential or no-slip), take the some pressure at the bottom (ideally the average), and prescribe this as P_0 in your equation. Does that work?
. I think this error message is incorrect, it is usually not the right thing to combine a free surface with velocity or traction boundary condition (if the free surface were at the top), but for your case it is the right thing to do. So I would just comment out the variable