Winkler foundation boundary condition?

[ricitron]

Is it possible to implement a Winkler foundation boundary condition for the bottom of a 2D box model? It would be a variable force on the bottom boundary equal to the boundary displacement times the gravity and mantle density. The idea is to simulate the flexure of the mechanical elastic lithosphere when a load is placed on the surface. Would this be possible with a the current boundary options or does it require a code modification?

[gassmoeller]

Hi @ricitron and sorry for the late reply. In theory our traction boundary condition could be changed to create a force that depends on other quantities (like the Winkler foundation). However, how difficult that will be depends on a few factors:

• our current interface for the boundary traction model is boundary_traction = boundary_traction(position, boundary_normal_vector) i.e. each boundary traction plugin receives the current location and normal vector and has to compute the resulting traction (=force/area)
• like you described when you want to implement the Winkler foundation you need a way to compute:
  • boundary displacement
  • gravity
  • density

For all of these there are 2 options:

1. you can compute (or approximate) the property directly from the input position, this is then easy to implement
2. you need to compute them accurately based on the current solution, this will be harder, because we are missing the correct interface (e.g. the material model gets a lot of information to compute its properties, but the boundary traction plugins do not)

Going with option 1 for a second: You could for example:

• compute boundary displacement as (original_boundary_location - position) if you know the original location of the boundary
• compute gravity from the gravity model, which uses the following interface: gravity = gravity(position) (since you know position in the boundary traction plugin everything works)
• approximate density with a known constant density

This would be relatively simple to implement and may even be possible without code modification using the function plugin. E.g. the following expression creates a traction in all directions (i.e. a variable pressure at the bottom boundary) that depends on the depth in a box model. Keep in mind you need to prescribe the full pressure, not just the dynamic pressure, and the function expression will still be multiplied with the normal of the boundary:

set Function constants = g = 10, rho = 3300, y0 = 100000
set Function expression = rho * g * (y0-y); rho * g * (y0-y)

If an approximate density is not enough, or you cannot compute the boundary displacement purely based on the current position of the boundary, then things get complicated. We would probably have to think about extending the interface, which is not impossible, but a bit of work.

I think it would be a cool extension to ASPECT to have an example or test of this. Let me know in case you want to contribute something like this.

[ricitron]

Thanks a ton for the reply. I will test this out (after AGU) and see how it goes with the option 1 you suggested.