Brittle thrust wedges benchmark
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. 203 and Fig. 204.
[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
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
# 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
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
set GMRES solver restart length = 200
# The discontinuous composition bound preserving limiter produces sharp
# interfaces between compositional layers.
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
# 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
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. 205). 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. 206). 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.