# Brittle thrust wedges benchmark#

*This section was contributed by Sibiao Liu, Stephanie Sparks, John Naliboff,
Cedric Thieulot, and Wolfgang Bangerth.*

Thrusting of brittle crust by applying compressive forces can lead to large deformations. The process is complicated to model because the rheology of cold, brittle crust is substantially more complicated than that of the hot, ductile rocks in the mantle. At the same time, the processes that act in such situations are surprisingly easy to replicate and visualize using simple “sand box” experiments in which one fills a volume with layers of differently-colored sand and compresses or stretches the volume. Examples of the patterns one can then observe in these do-it-yourself models are shown in Fig. 234 and Fig. 235.

[Buiter *et al.*, 2016] organized new comparison experiments
between these kinds of analogue and numerical models to investigate this kind
of brittle thrust wedge behavior. The benchmark here aims to verify that the
wedge models using follows other numerical results and the analytical wedge
theory shown in this paper. In particular, input files
(benchmarks/buiter_et_al_2016_jsg) are provided for reproducing the
numerical simulations of stable wedge experiment 1 and unstable wedge
experiment 2 with the same model setups.

A number of model sets of prescribed material behavior are required to simulate the brittle thrust formation. For example, although the material in the numerical model has a visco-plastic rheology, it performs plastic yielding at the beginning of shortening due to the non-viscous sand. We prescribe plastic strain-weakening behavior, with the internal angle of friction diminishing between total finite strain invariant values of 0.5 and 1.0, to mimic the softening from peak to dynamic stable strength which correlates with sand dilation.

In sandbox-type models, an important role is played by the boundaries and the frictional sliding of sand against these boundaries. For the top boundary condition, zero traction (“open”) and a sticky air layer is used to approximate a free surface. Additional testing revealed that using a true free surface leads to significant mesh distortion and associated numerical instabilities. We also apply a rigid block that approximates a mobile wall with a constant velocity of 2.5 cm/hour on the right-hand side boundary to drive the deformation in the sand layers. The following listing shows key portions of the parameter file that describes this kind of setup:

```
# Spatial domain of different compositional fields:
# Quartz sand (qsand) has two layers with 1 cm thickness each.
# Corundum sand (csand) is in the middle of the 2-layer quatz sand and is
# 1 cm thick.
# The boundary 2-mm-thin layers (bound) sit on both the bottom of the domain and
# between the sand and rigid indenter block. The right boundary has a constant
# inflow velocity, which pushes a rigid block that approximates a mobile wall.
# Movement of the rigid block drives deformation in the sand layers.
# The boundary between the rigid block and boundary layers produces a sharp
# velocity discontinuity that localizes brittle deformation.
# The sticky air layer is set on top of the sand layer and approximates a
# free surface.
# Velocity boundary conditions
subsection Boundary velocity model
set Zero velocity boundary indicators = bottom #no slip
set Tangential velocity boundary indicators = left #free slip
# right - material inflow with a velocity of 2.5 cm/hour at the height of
# the pushing block. The velocity linearly decreases from the base of the
# rigid block to 0 cm/hour at the base of the model.
set Prescribed velocity boundary indicators = right:function
subsection Function
set Variable names = x,y
set Function constants = cm=0.01, h=3600, th=0.002
set Function expression = if(y>th, -2.5*cm/h, -(y/th)*2.5*cm/h); 0
end
end
# The top boundary is open (zero traction), which allows the sticky air to
# flow freely through it as topography develops along the wedge. Additional
# testing revealed that using a true free surface leads to significant mesh
# distortion and associated numerical instabilities.
subsection Boundary traction model
set Prescribed traction boundary indicators = top: zero traction
end
```

Accurate solver convergence is always challenging to achieve in numerical
thrust wedge models with a high spatial resolution (ca. 1 mm node spacing) and
a large viscosity contrast. Here, we suggest that several parameters should be
considered carefully. First, the nonlinear and linear solver tolerances should
be sufficiently strict to avoid numerical instabilities. Second, we use the
discontinuous Galerkin method
(`set Use discontinuous composition discretization = true`

) to ensure that the
discontinuous composition bound preserving limiter produces sharp interfaces
between compositional layers. Lastly, we use the harmonic averaging scheme for
material and viscosity is required to achieve reasonable convergence behavior.
The relevant parameters are shown here:

```
# Note that the Linear/Nonlinear solver tolerance should be sufficiently strict
# to avoid numerical instabilities.
set Nonlinear solver scheme = single Advection, iterated Stokes
set Nonlinear solver tolerance = 1e-7
set Max nonlinear iterations = 100
subsection Solver parameters
subsection Stokes solver parameters
set Linear solver tolerance = 1e-8
set Number of cheap Stokes solver steps = 0
# A higher restart length makes the solver more robust for large viscosity
# contrasts.
set GMRES solver restart length = 200
end
end
# The discontinuous composition bound preserving limiter produces sharp
# interfaces between compositional layers.
subsection Discretization
set Use discontinuous composition discretization = true
subsection Stabilization parameters
set Use limiter for discontinuous composition solution = true
set Global composition maximum = 1, 1, 1, 1, 1, 100
set Global composition minimum = 0, 0, 0, 0, 0, 0
end
end
# Material properties
# Using harmonic material averaging is required to achieve reasonable
# convergence behavior, particulally when the viscosity contrast is large.
subsection Material model
set Material averaging = harmonic average
set Model name = visco plastic
subsection Visco Plastic
# The viscosity contrast is 10^7 here and any value higher than this causes
# divergence of the solver.
set Minimum viscosity = 1e5
set Maximum viscosity = 1e12
set Viscosity averaging scheme = harmonic
end
end
```

*Experiment 1* tests whether model wedges in the stable domain of critical
taper theory remain stable when translated horizontally. A quartz sand wedge
with a horizontal base and a surface slope of 20 degrees is pushed 4 cm
horizontally by inward movement of a mobile wall at the right boundary with a
velocity of 2.5 cm/hour (Fig. 236). The basal angle is zero
(horizontal), a thin layer separates the sand and boundary to ensure minimum
coupling between the wedge and bounding box during translation, and a sticky
air layer is used above the wedge. Further, the purely plastic material should
not undergo any deformation during translation.

*Experiment 2* tests how an unstable subcritical wedge deforms to reach the
critical taper solution. In this experiment, horizontal layers of sand undergo
10 cm shortening by inward movement of a mobile wall with a velocity of 2.5
cm/hour (Fig. 237). Model results show thrust wedge generation near
the mobile wall through a combination of mainly in-sequence forward and
backward thrusting. The strain field highlights several incipient shear zones
that do not always accumulate enough offset to become visible in the material
field. The pressure field of the model remains more or less lithostatic, with
lower pressure values in (incipient) shear zones.