Benchmark for subduction zone modeling

In early October 2002 a group of researchers met at the University of Michigan at Ann Arbor, MI to discuss the modeling of the thermal structure and dynamics of subduction zones At the workshop it became clear that the community could greatly benefit from a set of benchmarks that allow for code testing and comparisons. We can identify two fundamental approaches for subduction zone modeling: a) fully dynamic, where the deformation of the slab is computed using descriptions of rheology and buoyancy; and b) wedge dynamical models, where the geometry and velocity of the slab is imposed kinematically, with a dynamic solution only for the wedge. This draft paper formulates a set of benchmarks of increasing complexity focusing on the second category. The first category requires a fundamentally more difficult approach. For a discussion of a number of potential dynamic benchmarks see www.geobench.org. We hope that this this set of benchmarks will evolve to a standard in subduction modeling, similar to the role of the mantle convection benchmarks formulated in Blankenbach et al., 1989; Busse et al., 1993; and Van Keken et al., 1997.

Please contact Peter van Keken (keken at umich.edu; 1-734-764-1497) for further information.

Most recent description of benchmarks (December 2005)

Groups who participated in the benchmark development:
PGC: Pacific Geoscience Group: Claire Currie (now at Dalhousie), John He, Kelin Wang. Finite elements. claire.currie at dal.ca; wang at pgc.nrcan.gc.ca ; he at pgc.nrcan.gc.ca
UM: University of Michigan: Peter van Keken. Sepran finite element package. Linear: linear triangles with SUPG for heat equation, Taylor-Hood for Stokes equation; or Quadratic: quadratic triangles with SUPG for heat equation and penalty function method (Crouxiez-Raviart elements) for Stokes). keken at umich.edu
Purdue: Purdue University: Scott King. Conman. sking at purdue.edu
LDEO: Lamont Doherty: Rich Katz (now at Cambridge), Marc Spiegelman. Finite volume multigrid. mspieg at ldeo.columbia.edu; rfk22 at cam.ac.uk
Brown: Brown University: Amandine Cagnioncle, Jennifer Rilling, Marc Parmentier, Chad Hall. Finite elements for Stokes equation, Smolarkiewicz (J.Comp.Physics, 54, 1983) for heat equation. em_parmentier at brown.edu ; amandine_cagnioncle at brown.edu
Cardiff: Cardiff University: Huw Davies. Semi-analytical.
WHOI: Woods Hole: Mark Behn. Comsol 3.2b finite element code. mbehn at whoi.edu
NTU: National Taiwan University. Shu-Chuan Lin. Finite volume multigrid. skylin0 at ntu.edu.tw

Recent results (under construction; last update March 1, 2007)

We found it most relevant to use the temperature at the slab wedge interface to compare our models. We use spot temperatures at (54,54) and (210,210) and the L2 norm of the temperature along the slab (at 6 km intervals in depth) from 0 to 210. This allows us to plot the behavior of the model solution using a single parameters from the benchmark output.

Case 1a (updated 3/1/07)

Temperature at (54,54) as function of resolution near the corner point.
Temperature at (60,60) as function of resolution near the corner point.
Temperature at (210,210) as function of resolution near the corner point).
L2 norm of temperature from (0,0) to (210,210) as function of resolution near the corner point.
L2 norm of temperature in wedge spanned by (54,54), (120,54) and (120,120) as function of resolution near the corner point.
Filled symbols for UM represent results for quadratic elements for heat equation, compared to those with linear elements as open symbols.
Comparison by Amandine Cagnioncle for the Brown code with either 1 to 3 steps in the Smolarkiewicz scheme and a FD SUPG method (updated 9/17/06)
Temperature data files for
UM (highest resolution);
PGC (highest resolution);
PU (highest resolution);
LDEO: 0.5 km resolution 1 km. 2 km.
WHOI: 10 km global resolution, 0.5 km resolution in the nose; 1 km resolution in the nose; 5 km resolution in the nose;
NTU: 0.7 km resolution in the nose;
Brown: iord=1 0.4x0.4 grid; 1x1 grid; 2x2 grid; 3x3 grid;
Brown: iord=2 0.4x0.4 grid; 1x1 grid; 2x2 grid; 3x3 grid;
Brown: iord=3 0.4x0.4 grid; 1x1 grid; 2x2 grid; 3x3 grid;
In alternative format: Cardiff.

Case 1a_alt (updated 4/1/07)

Case with a 90 degree dipping slab with slab interface at x=150 km. All other parameters as in 1a. .

Temperature at (150,54) as function of resolution near the corner point.
Temperature at (150,60) as function of resolution near the corner point.
Temperature at (150,210) as function of resolution near the corner point.
L2 norm of temperature from (150,0) to (150,210) as function of resolution near the corner point.
L2 norm of temperature in square spanned by (150,54) and (216,120) as function of resolution near the corner point.
PDF with temperature plot for UM.
Temperature data files. Slab temperature at 54 km depth is about 400 (compared to 320 for 1a).
UM (300 m resolution in nose).
LDEO 0.5 km, LDEO 1.0 km, LDEO 2.0 km.

Case 1b (updated 4/25/07)

Temperature at (54,54) as function of resolution near the corner point.
Temperature at (60,60) as function of resolution near the corner point.
Temperature at (210,210) as function of resolution near the corner point.
L2 norm of temperature from (0,0) to (210,210) as function of resolution near the corner point.
L2 norm of temperature in wedge spanned by (54,54), (120,54) and (120,120) as function of resolution near the corner point.
The UM convergence results show good behavior and convergence toward the solution obtained in 1a for a moderate interval ('transition') in which the velocity boundary condition at the slab interface ramps up from zero to full velocity. Note the rather poor convergence behavior if such modification is not employed. To highlight the change compare the cases without this transition with the ones with transition. See also Appendix A in Van Keken et al., Geochem. Geophys. Geosyst., 2002.
'Best' temperature data files for UM; PGC; PU; LDEO. WHOI;

Dynamic pressure and pressure gradients in the wedge at two arbitrary points (point 1 corresponds to indices (12,10): x=66,z=54 and point 2 to indices (100,60): x=594,z=354) with pressure in MPa and pressure gradient in MPa/km:

                        P_1          P_2       dP/dx_1      dP/dy_1       dP/dx_2       dP/dy_2

analytical            -929.46       -7.34       55.97        12.71        0.02935       0.02748
PGC (best model)      -929.53       -7.35       51.16        11.45        0.03047       0.02758            
UM (linear; high res) -930.51       -7.70       56.24        12.51        0.0294        0.0275
UM (quad; high res)   -930.18       -7.35     
LDEO (0.5 km)         -916          -7.3        54.5         11.9         0.0293        0.0275
PU (360x180; 1.1km)   -927.8        -7.11     
NTU (0.88 km)         -911.9         2.2        55.63        16.02        0.02935       0.02751
WHOI                  -930.26       -7.35       56.86        14.45        0.0293        0.0277

Pressure plots of pressure (with P=0 at [660,600]); dP/dx; dP/dy.
Second invariant of strainrate tensor (UM, 150m resolution) in 1e-16/s data file; plot of log(strainrate).

Case 1b_alt (updated 3/14/07)

Case with a 90 degree dipping slab with slab interface at x=150 km. All other parameters as in 1b.

Case 1c (updated 3/1/07)

Temperature at (54,54) as function of resolution near the corner point.
Temperature at (60,60) as function of resolution near the corner point.
Temperature at (210,210) as function of resolution near the corner point.
L2 norm of temperature from (0,0) to (210,210) as function of resolution near the corner point.
L2 norm of temperature in wedge spanned by (54,54), (120,54) and (120,120) as function of resolution near the corner point.

The results are quite similar to those shown for 1b, suggesting that the change in velocity boundary conditions at in and outflow boundaries of the the wedge do not influence the dynamics and temperature field near the corner point significantly.
Temperature data files for UM; PGC; PU; LDEO. WHOI (10 km global resolution, 0.5 km resolution in the nose);
Second invariant of strainrate tensor (UM, 150m resolution) in 1e-16/s data file; plot of log(strainrate).

Case 2a (updated 4/1/07)

Benchmark for subduction zone modeling

Recent results (under construction; last update March 1, 2007)

Case 1a (updated 3/1/07)

Case 1a_alt (updated 4/1/07)

Case 1b (updated 4/25/07)

Case 1b_alt (updated 3/14/07)

Case 1c (updated 3/1/07)

Case 2a (updated 4/1/07)

Case 2a alt (updated 4/23/07)

Case 2b (updated 4/8/07)

AGU Fall 2005 poster (draft in PDF)

Previous and obsolete results