Van Keken et al. (2008) wedge flow benchmark problem

[Reverie]

Hi all,
I was following the wedge flow benchmark cookbook and found that it is not yet finished. Running the provided prm file with the custom mesh gives a subducting slab much cooler than the benchmark:

In the right panel above the temperature contour above the slab interface is too sparse compared to the benchmark (around 390℃ or 663 K, see van Keken et al., 2008), the temperature at (60, 60km) on the slab is only 441.37 K (after resampling to a 111*101 regular mesh) :

I then returned to the regular rectangular mesh, and contrary to what is said in the cookbook, the Stokes equation converges and reaches a steady state (perhaps it was because I was using deal.II-v9.7.0, and on another device with deal.II-v9.6.1 it failed to converge with the same prm file):

The temperature contour near the slab is more compact now, but still not enough, T(60, 60km) = 503.90 K.

Then I employed a velocity ramp near the tip, as suggested in the benchmark paper, but this only raises the temperature at that point for several tens of degrees, still a long way to go. The only effective approach seems to be increasing the global mesh refinement…when the mesh spacing is ~1km, T(60, 60km) = 507.9K and when the spacing is further reduced (~500m) it became 550.99K, a mesh spacing of ~250m gives ~570K.

I also calculated the slab temperature and wedge temperature for the ~500m mesh spacing case, T_slab=681.73K and T_wedge=1056.55 K, both lower than the benchmark values (by about 100 degrees, similar to the offset of T(60, 60km))

The convergence with increasing refinement is expected, but the convergence rate is impractically slow. I tried to use AMR as well but cannot get the Stokes equation converge at the first timestep. The velocity discrepancy with respect to the analytical solution is like (for the ~500m mesh spacing case, and I imposed analytical velocity at the right boundary):

So the main issue is still the velocity field discrepancy near the wedge tip, even though the area is not that large now. I don’t know what to do now, any suggestions on how to further improve the velocity field (so temperature field as well) in this region would be very helpful.

Thanks in advance!

Here is the parameter file I used to calculate the ~500m mesh spacing case.
original.prm (7.0 KB)

[daniel.douglas]

Hi @Reverie

Thank you for taking the time to formally compare this cookbook to the actual benchmark values. As you point out this was something that I had meant to do as a follow up to the initial implementation of this cookbook, but I never got around to it. It’s also great that it converges with the default ASPECT mesh.

Besides increasing the resolution of the mesh, which you’ve already done, the other parameter I would try changing is the solver tolerances and/or the discretization. Maybe a stricter tolerance will help chip away at the error?

After looking at the .prm file you link and comparing it against van keken et al., 2008, one other thing you might try is making the thermal expansion coefficient 0. They seem to be using a constant density, and if I’m not mistaken the variable density in the models you’ve been running will lead to a variable thermal diffusivity which may be contributing to the temperature error as well. Thanks again for tackling this, looking forward to hearing if any of this helps!

Cheers,

Daniel

[maxrudolph]

Hi Daniel and @Reverie,

I continue to be interested in the problem of implementing corner flow models in ASPECT. In case you have not seen it, there is a very old pull request here:

[WIP] Vankeken et al. 2008 subduction benchmark

My pull request built off of a prescribed velocity plugin written by Jonathan Perry-Houts and Daniel built on what I had done. I could never get the velocity and temperature fields to agree with what is in the van keken et al. paper. Solver convergence was not a problem and if you suspect that it is, you could likely use a direct solver for a 2D problem like this. When I was trying to reproduce the benchmark, Peter van Keken sent me the benchmark velocity and temperature values on a regularly sampled grid. If nothing else, this might help you to isolate the discrepancy between ASPECT and the other codes.

Max

[Reverie]

Hi Daniel and Max @maxrudolph,

Thanks for your suggestions. Sorry for taking so long to reply. I tried them but there was no luck. Below are the T(60, 60km) values in different cases, all have ~500m mesh spacing:

Reference: 550.99K
Zero thermal expansivity: 549.776K;
Lower tolerance (Temperature solver tolerance = 1e-14, also zero thermal expansivity): 549.776K.

So the effects of thermal expansivity and solver tolerance seem to be minor here. Increasing the temperature polynomial degree by 1 (so Q3 was used for the temperature field) didn’t affect the result much.

I am also wondering if a finer initial mesh + lower degree of global refinement can be different from a coarse initial mesh + high degree of global refinement (like what I am doing here, a 11*10 initial mesh and global refinement = 5/6 etc.), and have just submitted runs with finer initial mesh and lower CFL.

@maxrudolph Max, thanks for your link to the PR almost 9 years ago! From your results of case1a there, the entire slab is colder than van Keken’s, but what is the unit of the difference there? Does that mean ASPECT cannot easily reproduce the temperature field for this problem even if the entire velocity field is prescribed? If the unit is Kelvin, then my case1b recreation seems far worse (around 100K descrepancy), I will take a closer look at your prm file later.

I was planning to incorporate more physics based on this wedge flow model, so that is why I am trying to reproduce the benchmark now. Maybe I should just use real-world rheology/geometry to bypass this mathematical singularity.

Best,
Hang