diff --git a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_t.prm b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_t.prm
index ca7b64af468..97b4efe6cb4 100644
--- a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_t.prm
+++ b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_t.prm
@@ -19,6 +19,7 @@ set Additional shared libraries = libsimple_nonlinear.so
 subsection Solver parameters
   subsection Stokes solver parameters
     set Linear solver tolerance = 1e-5
+    set Stokes solver type =  block AMG
   end
 
   subsection Newton solver parameters
diff --git a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_v.prm b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_v.prm
index 1ceef140ad8..424ec03d84d 100644
--- a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_v.prm
+++ b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_v.prm
@@ -19,6 +19,7 @@ set Additional shared libraries = libsimple_nonlinear.so
 subsection Solver parameters
   subsection Stokes solver parameters
     set Linear solver tolerance = 1e-15
+    set Stokes solver type = block AMG
   end
 
   subsection Newton solver parameters
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/README.md b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/README.md
index 95b5d298295..b82929c0935 100644
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/README.md
+++ b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/README.md
@@ -1,11 +1,19 @@
 # Newton Solver Benchmark Set - Spiegelman et at. (2016)
 
 The files in [this directory](https://github.com/geodynamics/aspect/tree/main/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016)
-can be used to recreate the Spiegelman et al. (2016) flow figures from
-{cite:t}`fraters:etal:2019`. Adjust both the `metabash.sh` and `bash.sh` files to
-reflect the parameter search you are interested in and run the metabash script.
-If you have run exactly the same runs as shown in the paper, you can use the
-`cp_P-minLT-res-it-vel_BV.sh` script to recreate the figures in that paper.
-This script requires that the drucker prager compositions plugin is installed.
-It should be possible to replicate the same behavior with the visco plastic
-plugin.
+can be used to recreate the figures from
+{cite:t}`fraters:etal:2019` reproducing the {cite:t}`spiegelman:etal:2016` benchmark.
+`input.prm` is the setup of one of the benchmark
+models and `run.sh` will execute a series of model runs. Before you can execute
+the benchmark you will have to compile the plugin `drucker_prager_compositions.cc` in
+the same directory, the process of compiling plugins is described in [](sec:benchmark-run).
+After you have executed `run.sh`, you can use the gnuplot plotting script
+`plot.gnuplot` to recreate Fig. 4 of {cite:t}`fraters:etal:2019`.
+
+```{figure-md} fig:benchmark-newton-spiegelman-2016
+<img src="doc/figure_4.png" />
+
+Convergence history of several models of the Spiegelman et al. benchmark: a reproduction of three of the nine pressure dependent Drucker–Prager cases with a resolution of 64 × 16 cells. Top: results for computations where linear systems are solved with a maximum relative tolerance of 0.9. Bottom: The same models solved with a tolerance of $10^{-8}$. The initial Picard iteration is always solved to a linear tolerance of $10^{-14}$. Left to right: different prescribed velocities of u0 = 2.5, 5, and 12.5 cm\yr, and and different background reference viscosities of respectively $\eta_{ref}$ = $10^{23}$, $10^{24}$ and $10^{25}$ Pa s. Horizontal axis: number of the non-linear (outer) iteration; and vertical axis: non-linear residual. "DC Picard" refers to a defect correction Picard iteration, see the paper describing this solver.
+```
+
+The nonlinear convergence behavior of ASPECT's different nonlinear solvers is shown in {numref}`fig:benchmark-newton-spiegelman-2016`.
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/auto_logs/.gitignore b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/auto_logs/.gitignore
deleted file mode 100644
index 5e7d2734cfc..00000000000
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/auto_logs/.gitignore
+++ /dev/null
@@ -1,4 +0,0 @@
-# Ignore everything in this directory
-*
-# Except this file
-!.gitignore
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/bash.sh b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/bash.sh
deleted file mode 100755
index 2b4603255fc..00000000000
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/bash.sh
+++ /dev/null
@@ -1,162 +0,0 @@
-#!/bin/bash
-version="test"
-build_dir="../../../build/"
-declare -a grid=(5) 
-declare -a UFS=("false" "true") # use failsafe
-declare -a NSP=("SPD") # "Symmetric" "PD" "none") #Use Newton stabilization preconditioner
-declare -a NSA=("SPD" "symmetric" "none") # "Symmetric" "PD" "none") #Use Newton stabilization A block
-declare -a agrid=(0)
-declare -a n=(50) #1 2 5 10 25 50 100)
-I=150
-declare -a P=(0 5 10 50 15 150 20 25 30) #0 5 10 15 20 150) # 1 2 3 4 5 10 15 25 50 100)
-declare -a LS=(50)
-declare -a ST=(1e-20)
-declare -a LT=("1e-5") # "1e-8") # "1e-8") #"1e-5") # "1e-8")
-declare -a NLT=("1e-14") # "1e-9" "1e-10" "1e-14" "1e-20")
-declare -a ABT=("1e-2")
-declare -a RSM=("true") #"true" #"false"
-declare -a SF=("9e-1")
-declare -a UDS=("true")
-declare -a AEW=("false") # always use Eisenstat Walker, even for picard
-AV=-1
-COMP="0" #"4e-12"
-declare -a OS=("9e-1" "1e-1" "1e-2" "1e-8") # "1e-3" "1e-4" "1e-5" "1e-6" "1e-6" "1e-7" "1e-8")
-materialmodelnameShort="DP" #"VP2"
-materialmodelname="DP"
-phi=30 # 0
-declare -a vel=(25 50 125) # 50 125) #25 50 125)
-declare -a BV=("1e23" "1e24" "5e24") #"1e23" "1e24" "5e24")
-processes=4
-
-CohesionLine="s/set Cohesions.*/      set Cohesions                              = 1e8,0/g"
-PhiLine="s/set List of angles of internal friction.*/      set List of angles of internal friction            = $phi,0/g"
-if [ $materialmodelnameShort == "DP" ]; then
- materialmodelname="drucker prager compositions"
-elif [ $materialmodelnameShort == "vM" ]; then
- materialmodelname="drucker prager compositions"
-elif [ $materialmodelnameShort == "SNL" ]; then
- materialmodelname="simple nonlinear compositions"
- CohesionLine="s/ set List of cohesion of fields.*/ /g"
- PhiLine="s/set List of angles of internal friction of fields.*/ /g"
-elif [ $materialmodelnameShort == "VP2" ]; then
- materialmodelname="viscoplastic2"
-fi
-
-SOLVER_SHORT="NS" #"itAdandSt" #"NS" #NS"
-SOLVER="NS"
-if [ $SOLVER_SHORT == "NS" ]; then
- SOLVER="iterated Advection and Newton Stokes"
-elif [ $SOLVER_SHORT == "itAdandSt" ]; then
- SOLVER="iterated Advection and Stokes"
-fi
-
-
-for i_grid in "${grid[@]}"
-do
- for i_agrid in "${agrid[@]}"
- do
-  for i_NSP in "${NSP[@]}"
-  do
-   for i_NSA in "${NSA[@]}"
-   do
-    for i_NLT in "${NLT[@]}"
-    do
-     for i_ABT in "${ABT[@]}"
-     do
-      for i_ST in "${ST[@]}"
-      do
-       for i_UDS in "${UDS[@]}"
-       do
-        for i_SF in "${SF[@]}"
-        do
-         for i_LT in "${LT[@]}"
-         do
-          for i_n in "${n[@]}"
-          do
-           for i_UFS in "${UFS[@]}"
-           do
-            for i_P in "${P[@]}"
-            do
-             for i_LS in "${LS[@]}"
-             do
-              for i_AEW in "${AEW[@]}"
-              do
-               for i_OS in "${OS[@]}"
-               do
-                for i_RSM in "${RSM[@]}"
-                do
-                 for i_BV in "${BV[@]}"
-                 do
-                  for i_vel in "${vel[@]}"
-                  do
-                  U="7.92219116e-11"
-                  if [ $i_vel == 25 ]; then
-                   U="7.92219116e-11"
-                  elif [ $i_vel == 50 ]; then
-                   U="1.58443823e-10"
-                  elif [ $i_vel == 125 ]; then
-                   U="3.96109558e-10"
-                  fi
-                  
-                  dirname_clean="$version""$materialmodelnameShort""_""$SOLVER_SHORT""_ST""$i_ST""_UFS""$i_UFS""_NSP-""$i_NSP""_NSA-""$i_NSA""_C""$COMP""_g""$i_grid""_ag""$i_agrid""_AEW""$i_AEW""_UDS""$i_UDS""_SF""$i_SF""_NLT""$i_NLT""_ABT""$i_ABT""_LT""$i_LT""_mLT""$i_OS""_I""$I""_P""$i_P""_EW1""_theta1""_LS""$i_LS""_RSM""$i_RSM""_AV""$AV""_phi""$phi""_vel""$i_vel""_BV""$i_BV""_n""$i_n"
-                  dirname="results/""$dirname_clean"
-                  infilename="$dirname""/input.prm"
-                  outfilename="$dirname""/output.log"
-                  errorfilename="$dirname""/error.log"
-                  outplotfilename="$dirname""/plot.dat"
-                  
-                  mkdir -p $dirname
-                  
-                  echo "$dirname"
-                  sed  \
-                   -e "$PhiLine" \
-                   -e "$CohesionLine" \
-                   -e "s/set Function expression = if(x<60e3,.*/    set Function expression = if(x<60e3,$U,-$U);0/g" \
-                   -e "s/set Nonlinear Newton solver switch tolerance.*/     set Nonlinear Newton solver switch tolerance = $i_ST/g" \
-                   -e "s/set Reference viscosity .*/    set Reference viscosity = $i_BV/g" \
-                   -e "s/set Output directory .*/set Output directory = results\/$dirname_clean/g" \
-                   -e "s/set Model name .*/    set Model name = $materialmodelname/g" \
-                   -e "s/set Nonlinear solver scheme.*/set Nonlinear solver scheme = $SOLVER/g" \
-                   -e "s/set Initial global refinement .*/  set Initial global refinement          = $i_grid/g" \
-                   -e "s/set Initial adaptive refinement .*/  set Initial adaptive refinement        = $i_agrid/g" \
-                   -e "s/set Stress exponents of fields .*/      set List of stress exponents of fields                       = $i_n, 1/g" \
-                   -e "s/set Viscosity averaging p .*/    set Viscosity averaging p = $AV/g" \
-                   -e "s/set Max nonlinear iterations .*/set Max nonlinear iterations = $I /g" \
-                   -e "s/set Linear solver tolerance .*/set Linear solver tolerance = $i_LT/g" \
-                   -e "s/set Nonlinear solver tolerance.*/set Nonlinear solver tolerance = $i_NLT/g" \
-                   -e "s/set Linear solver A block tolerance.*/set Linear solver A block tolerance = $i_ABT/g" \
-                   -e "s/set Reference compressibility .*/    set Reference compressibility = $COMP/g" \
-                   -e "s/set Max pre-Newton nonlinear iterations .*/    set Max pre-Newton nonlinear iterations = $i_P/g" \
-                   -e "s/set Use Newton failsafe .*/set Use Newton failsafe = $i_UFS/g" \
-                   -e "s/set Stabilization preconditioner .*/set Stabilization preconditioner = $i_NSP/g" \
-                   -e "s/set Stabilization velocity block .*/set Stabilization velocity block = $i_NSA/g" \
-                   -e "s/set SPD safety factor .*/set SPD safety factor = $i_SF/g" \
-                   -e "s/set Use deviator of strain-rate .*/set Use deviator of strain-rate = $i_UDS/g" \
-                   -e "s/set Max Newton line search iterations .*/    set Max Newton line search iterations = $i_LS/g" \
-                   -e "s/set Maximum linear Stokes solver tolerance .*/set Maximum linear Stokes solver tolerance = $i_OS/g" \
-                   -e "s/set Use Newton residual scaling method .*/    set Use Newton residual scaling method = $i_RSM/g" \
-                  input.prm > "$infilename"
-                  
-                  nohup mpirun -np $processes $build_dir./aspect $infilename > $outfilename 2>$errorfilename
-                  
-                  grep "Relative nonlinear residual" $outfilename > $outplotfilename
-                  
-                  done
-                 done
-                done
-               done
-              done
-             done
-            done
-           done
-          done
-         done
-        done
-       done
-      done
-     done
-    done
-   done
-  done
- done
-done
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/bash_autos/.gitignore b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/bash_autos/.gitignore
deleted file mode 100644
index 5e7d2734cfc..00000000000
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/bash_autos/.gitignore
+++ /dev/null
@@ -1,4 +0,0 @@
-# Ignore everything in this directory
-*
-# Except this file
-!.gitignore
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/cp_P-minLT-res-it-vel_BV.sh b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/cp_P-minLT-res-it-vel_BV.sh
deleted file mode 100755
index 9c7a59352fa..00000000000
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/cp_P-minLT-res-it-vel_BV.sh
+++ /dev/null
@@ -1,258 +0,0 @@
-#!/bin/bash
-declare -a version="c3"
-#declare -a boundary_type=("velocities") # "tractions")
-#for i_boundary_type in "${boundary_type[@]}"
-#do
-input_file="input_""$i_boundary_type"".prm"
-declare -a grid=(4) #4 5 6 7 #2 3 4 5)
-agrid=0
-declare -a NSP=("SPD") # "symmetric" "PD" "none") #Use Newton stabilization preconditioner
-declare -a NSA=("SPD") # "symmetric" "none") # "Symmetric" "PD" "none") #Use Newton stabilization A block
-declare -a UFS=("false" "true")
-declare -a n=(50) # 5) #3 5 #1 2 5 10 25 50 100)
-I=150
-declare -a P=(0 5 10 15 20 25 30 50 150) # 5 10 100) #0 5 10 100 
-declare -a ST=(1e-20)
-declare -a LS=(0 5 10) # 5 50) # 5 10 50) # 5 10) # 0,5,10
-declare -a LT=("1e-5") # "1e-8") # 6 7 8) # 5,6,7,8
-declare -a NLT=("1e-14")
-declare -a RSM=("false" "true")
-declare -a UDS=("true") #"true" #"false"
-declare -a SF=("9e-1" "9.9e-1") # "9.99e-1" "1")
-declare -a AEW=("false")
-AV=-1
-ATB="1e-2"
-COMP="0" #"4e-12"
-materialmodelnameShort="DP" #"VP2"
-materialmodelname="DP"
-phi=30
-declare -a vel=(25 50 125)
-declare -a BV=("1e23" "1e24" "5e24")
-SOLVER_SHORT="NS" #"itAdandSt"
-SOLVER="iterated Advection and Stokes"
-if [ $SOLVER_SHORT == "NS" ]; then
-SOLVER="iterated Advection and Newton Stokes"
-elif [ $SOLVER_SHORT == "itAdandSt" ]; then
-SOLVER="iterated Advection and Stokes"
-fi  
-declare -a OS=("9e-1" "1e-1" "1e-2" "1e-8") #"1e-7" "1e-6" "1e-5" "1e-4" "1e-3" "1e-2" "1e-1" "5e-1" "9e-1") #"9e-1" "5e-1" "1e-1" "1e-2" "1e-4")
-
-for i_UFS in "${UFS[@]}"
-do
-for i_NSP in "${NSP[@]}"
-do
-for i_NSA in "${NSA[@]}"
-do
-for i_ST in "${ST[@]}"
-do
-for i_RSM in "${RSM[@]}"
-do
-for i_grid in "${grid[@]}"
-do
-for i_LT in "${LT[@]}"
-do
-for i_LS in "${LS[@]}"
-do
-for i_UDS in "${UDS[@]}"
-do
-for i_SF in "${SF[@]}"
-do
-for i_AEW in "${AEW[@]}"
-do
-for i_NLT in "${NLT[@]}"
-do
-for i_n in "${n[@]}"
-do
-
-basename="$version""P-minLT-res-it-vel-BV_g""$i_grid""_ag""$agrid""_ST""$i_ST""_UFS""$i_UFS""_NSP-""$i_NSP""_NSA-""$i_NSA""_ATB""$ATB""_n""$i_n""_I""$I""_LS""$i_LS""_LT""$i_LT""_NLT""$i_NLT""_RSM""$i_RSM""_""$SOLVER_SHORT""_AV""$AV""_phi""$phi""_comp""$COMP""_UDS""$i_UDS""_SF""$i_SF""_AEW""$i_AEW"
-#final_plot_file_name=${final_file_name::-4}".gnuplot"
-final_plot_file_name="image_plotfiles/""$basename"".gnuplot"
-final_png_file_name="image_collections/$basename.png"
-rm $final_plot_file_name
-
-echo "set terminal pngcairo size 3000,2400 enhanced font 'Verdana,20'" >> $final_plot_file_name
-echo "set output '$final_png_file_name'" >> $final_plot_file_name
-echo "set multiplot title \"{/*2.0 Final nonlinear iteration and final nonlinear residual as function of grid resolution and Max nonlinar tolerance}\\n\\n{/*1.5 Grid: $i_grid, LS: $i_LS, ST: $i_ST, NLT: $i_NLT, RSM: $i_RSM, phi: $phi, comp: $COMP, solver: $SOLVER_SHORT, UFS: $i_UFS, NSP: $i_NSP, NSA: $i_NSA, UDS: $i_UDS, SF: $i_SF, AEW: $i_AEW}\"" >> $final_plot_file_name
-echo "set rmargin 2" >> $final_plot_file_name
-echo "set bmargin 2.2" >> $final_plot_file_name
-for i_BV in "${BV[@]}"
-do
-for i_vel in "${vel[@]}"
-do
-
-final_file_name="image_data/""$basename""_vel""$i_vel""_BV""$i_BV"".dat" #P-NLT-res-it_g""$i_grid""_n""$i_n""_I""$I""_LS""$i_LS""_LT""$i_LT""_minLT""$i_OS""_""$SOLVER_SHORT"".dat"
-rm $final_file_name
-
-
-for i_P in "${P[@]}"
-do
-for i_OS in "${OS[@]}"
-do
-#dirname="b3_""$i_boundary_type""_g""$i_grid""_n""$i_n""_I""$I""_P""$i_P""_LS""$i_LS""_LT""$i_LT""_NLT""$i_NLT""_UA""$i_UA""_minLT""$i_OS""_""$SOLVER_SHORT""_phi""$phi""_vel""$i_vel""_BV""i_BV"
-dirname="results/""$version""$materialmodelnameShort""_""$SOLVER_SHORT""_ST""$i_ST""_UFS""$i_UFS""_NSP-""$i_NSP""_NSA-""$i_NSA""_C""$COMP""_g""$i_grid""_ag""$agrid""_AEW""$i_AEW""_UDS""$i_UDS""_SF""$i_SF""_NLT""$i_NLT""_ABT""$ATB""_LT""$i_LT""_mLT""$i_OS""_I""$I""_P""$i_P""_EW1""_theta1""_LS""$i_LS""_RSM""$i_RSM""_AV""$AV""_phi""$phi""_vel""$i_vel""_BV""$i_BV""_n""$i_n"
-
-infilename="$dirname""/input.prm"
-outfilename="$dirname""/output.log"
-outplotfilename="$dirname""/plot.dat"
-#infilename="g""$i_grid""_n""$i_n""_""$I""I_""$i_P""P_LS""$LS""_LT""$LT""/input.prm"
-#outfilename="g""$i_grid""_n""$i_n""_""$I""I_""$i_P""P_LS""$LS""_LT""$LT""/output.log"
-#outplotfilename="g""$i_grid""_n""$i_n""_""$I""I_""$i_P""P_LS""$LS""_LT""$LT""/plot.dat"
-#dirname="g""${grid[i_grid]}""_n""${n[i_n]}""_""$I""I_""${P[i_P]}""P_LS""$LS""_LT""$LT"
-
-#mkdir -p $dirname
-
-echo "$dirname"
-
-awk -v dirname="$dirname" -v var1="$i_P" -v var2="$i_OS" 'END{print dirname ", " var1 ", " var2 ", " substr($11,1,length($11)-1) ", " substr($10,1,length($10)-1)}' $outplotfilename >> $final_file_name
-#sed  \
-#-e "s/set Nonlinear solver scheme.*/set Nonlinear solver scheme = $SOLVER/g" \
-#$input_file > "$infilename"
-
-#mpirun -np 1 ../../.././aspect $infilename > $outfilename
-#../../.././aspect $infilename > $outfilename
-
-#grep "Total relative residual after nonlinear iteration " $outfilename > $outplotfilename
-done
-done
-if [ "$i_vel" == "25" ] && [ "$i_BV" == "1e23" ]; then
-   echo "set origin 0,0.64" >> $final_plot_file_name
-   echo "set colorbox vert user origin .93,.0275 size .02,.92" >> $final_plot_file_name
-   echo "set xlabel \"Maximum Picard iterations, vel = 25, BV = 1e23\" offset 0,0.75" >> $final_plot_file_name
-elif [ "$i_vel" == "50" ] && [ "$i_BV" == "1e23" ]; then
-   echo "set origin 0.31,0.64" >> $final_plot_file_name
-   echo "unset colorbox" >> $final_plot_file_name
-   echo "set xlabel \"Maximum Picard iterations, vel = 50, BV = 1e23\" offset 0,0.75" >> $final_plot_file_name
-elif [ "$i_vel" == "125" ] && [ "$i_BV" == "1e23" ]; then
-   echo "set origin 0.62,0.64" >> $final_plot_file_name
-   echo "unset colorbox" >> $final_plot_file_name
-   echo "set xlabel \"Maximum Picard iterations, vel = 125, BV = 1e23\" offset 0,0.75" >> $final_plot_file_name
-elif [ "$i_vel" == "25" ] && [ "$i_BV" == "1e24" ]; then
-   echo "set origin 0,0.32" >> $final_plot_file_name
-   echo "unset colorbox" >> $final_plot_file_name
-   echo "set xlabel \"Maximum Picard iterations, vel = 25, BV = 1e24\" offset 0,0.75" >> $final_plot_file_name
-elif [ "$i_vel" == "50" ] && [ "$i_BV" == "1e24" ]; then
-   echo "set origin 0.31,0.32" >> $final_plot_file_name
-   echo "unset colorbox" >> $final_plot_file_name
-   echo "set xlabel \"Maximum Picard iterations, vel = 50, BV = 1e24\" offset 0,0.75" >> $final_plot_file_name
-elif [ "$i_vel" == "125" ] && [ "$i_BV" == "1e24" ]; then
-   echo "set origin 0.62,0.32" >> $final_plot_file_name
-   echo "unset colorbox" >> $final_plot_file_name
-   echo "set xlabel \"Maximum Picard iterations, vel = 125, BV = 1e24\" offset 0,0.75" >> $final_plot_file_name
-elif [ "$i_vel" == "25" ] && [ "$i_BV" == "5e24" ]; then
-   echo "set origin 0,0" >> $final_plot_file_name
-   echo "unset colorbox" >> $final_plot_file_name
-   echo "set xlabel \"Maximum Picard iterations, vel = 25, BV = 5e24\" offset 0,0.75" >> $final_plot_file_name
-elif [ "$i_vel" == "50" ] && [ "$i_BV" == "5e24" ]; then
-   echo "set origin 0.31,0" >> $final_plot_file_name
-   echo "unset colorbox" >> $final_plot_file_name
-   echo "set xlabel \"Maximum Picard iterations, vel = 50, BV = 5e24\" offset 0,0.75" >> $final_plot_file_name
-elif [ "$i_vel" == "125" ] && [ "$i_BV" == "5e24" ]; then
-   echo "set origin 0.62,0" >> $final_plot_file_name
-   echo "unset colorbox" >> $final_plot_file_name
-   echo "set xlabel \"Maximum Picard iterations, vel = 125, BV = 5e24\" offset 0,0.75" >> $final_plot_file_name
-fi
-echo "set size 0.29,0.32" >> $final_plot_file_name
-echo "" >> $final_plot_file_name
-echo "set ytics autofreq" >> $final_plot_file_name
-echo "set ylabel \"Maximum nonlinear tolerance\" offset 1,0" >> $final_plot_file_name
-echo "set cblabel \"Final nonlinear residual (square)\"" >> $final_plot_file_name
-echo "set palette rgb 33,13,10" >> $final_plot_file_name
-echo "set logscale y" >> $final_plot_file_name
-echo "set logscale cb" >> $final_plot_file_name
-echo "set format cb '%.0e'" >> $final_plot_file_name
-echo "set cbtics nomirror" >> $final_plot_file_name
-echo "" >> $final_plot_file_name
-echo "set format y '%.0e'" >> $final_plot_file_name
-echo "set xtics 10" >> $final_plot_file_name
-echo "set mxtics 2" >> $final_plot_file_name
-echo "set xrange [-3:53]" >> $final_plot_file_name
-echo "set yrange [1e-9:10]" >> $final_plot_file_name
-echo "set cbrange [1e-20:1e-1]" >> $final_plot_file_name
-echo "unset key" >> $final_plot_file_name
-echo "plot \"""$final_file_name""\" u 2:3:(3):4 w p ls 5 ps 10 lc palette " >> $final_plot_file_name
-#echo "set colorbox horiz user origin .1,.07 size .7,.04" >> $final_plot_file_name
-#echo "set cblabel \"Final nonlinear iteration (circle)\" offset 0,0.75" >> $final_plot_file_name
-echo "unset colorbox" >> $final_plot_file_name
-echo "unset logscale cb" >> $final_plot_file_name
-echo "set format cb '%.0f'" >> $final_plot_file_name
-echo "set cbtics nomirror" >> $final_plot_file_name
-echo "set style fill solid" >> $final_plot_file_name
-echo "set palette rgb 33,13,10" >> $final_plot_file_name
-echo "set cbrange [0:100]" >> $final_plot_file_name
-echo "set cbtics 20" >> $final_plot_file_name
-echo "set mcbtics 4" >> $final_plot_file_name
-#echo "plot \"""$final_file_name""\" u 2:3:(3):5 w p ls 7 ps 6 lc palette " >> $final_plot_file_name
-#echo "plot \"""$final_file_name""\" u 2:3:(3):5 w p ls 6 ps 6 lc rgb \"black\" " >> $final_plot_file_name
-echo "plot \"""$final_file_name""\" u 2:3:(sprintf('%.0f',\$5)) w labels font 'Times Bold,24' " >> $final_plot_file_name
-#echo "plot \"""$final_file_name""\" u 2:3:(2) w circle lc rgb \"black\" " >> $final_plot_file_name
-#echo "plot \"""$final_file_name""\" u 2:3:(1.5):5 w circle lc palette " >> $final_plot_file_name
-#echo "pause -1" >> $final_plot_file_name
-
-if [ "$i_vel" == "25" ] && [ "$i_BV" == "1e23" ]; then
-   echo "set origin 0.265,0.64" >> $final_plot_file_name
-elif [ "$i_vel" == "50" ] && [ "$i_BV" == "1e23" ]; then
-   echo "set origin 0.575,0.64" >> $final_plot_file_name
-elif [ "$i_vel" == "125" ] && [ "$i_BV" == "1e23" ]; then
-   echo "set origin 0.885,0.64" >> $final_plot_file_name
-   echo "unset colorbox" >> $final_plot_file_name
-elif [ "$i_vel" == "25" ] && [ "$i_BV" == "1e24" ]; then
-   echo "set origin 0.265,0.32" >> $final_plot_file_name
-elif [ "$i_vel" == "50" ] && [ "$i_BV" == "1e24" ]; then
-   echo "set origin 0.575,0.32" >> $final_plot_file_name
-elif [ "$i_vel" == "125" ] && [ "$i_BV" == "1e24" ]; then
-   echo "set origin 0.885,0.32" >> $final_plot_file_name
-elif [ "$i_vel" == "25" ] && [ "$i_BV" == "5e24" ]; then
-   echo "set origin 0.265,0" >> $final_plot_file_name
-elif [ "$i_vel" == "50" ] && [ "$i_BV" == "5e24" ]; then
-   echo "set origin 0.575,0" >> $final_plot_file_name
-elif [ "$i_vel" == "125" ] && [ "$i_BV" == "5e24" ]; then
-   echo "set origin 0.885,0" >> $final_plot_file_name
-fi
-
-echo "set size 0.05,0.32" >> $final_plot_file_name
-echo "" >> $final_plot_file_name
-echo "unset ylabel" >> $final_plot_file_name
-echo "unset xlabel" >> $final_plot_file_name
-echo "set cblabel \"Final nonlinear residual (square)\"" >> $final_plot_file_name
-echo "set palette rgb 33,13,10" >> $final_plot_file_name
-echo "set logscale y" >> $final_plot_file_name
-echo "set logscale cb" >> $final_plot_file_name
-echo "set format cb '%.0e'" >> $final_plot_file_name
-echo "set cbtics nomirror" >> $final_plot_file_name
-echo "" >> $final_plot_file_name
-echo "set format y '%.0e'" >> $final_plot_file_name
-echo "set xtics 10" >> $final_plot_file_name
-echo "set mxtics 2" >> $final_plot_file_name
-echo "unset ytics" >> $final_plot_file_name
-echo "set xrange [147:153]" >> $final_plot_file_name
-echo "set yrange [1e-9:10]" >> $final_plot_file_name
-echo "set cbrange [1e-20:1e-1]" >> $final_plot_file_name
-echo "unset key" >> $final_plot_file_name
-echo "plot \"""$final_file_name""\" u 2:3:(3):4 w p ls 5 ps 10 lc palette " >> $final_plot_file_name
-echo "unset colorbox" >> $final_plot_file_name
-echo "unset logscale cb" >> $final_plot_file_name
-echo "set format cb '%.0f'" >> $final_plot_file_name
-echo "set cbtics nomirror" >> $final_plot_file_name
-echo "set style fill solid" >> $final_plot_file_name
-echo "set palette rgb 33,13,10" >> $final_plot_file_name
-echo "set cbrange [0:100]" >> $final_plot_file_name
-echo "set cbtics 20" >> $final_plot_file_name
-echo "set mcbtics 4" >> $final_plot_file_name
-echo "plot \"""$final_file_name""\" u 2:3:(sprintf('%.0f',\$5)) w labels font 'Times Bold,24' " >> $final_plot_file_name
-
-gnuplot $final_plot_file_name
-done
-done
-done
-done
-done
-done
-done
-done
-done
-done
-done
-done
-done
-done
-done
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/doc/figure_4.png b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/doc/figure_4.png
new file mode 100644
index 00000000000..79a2a7f7b54
Binary files /dev/null and b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/doc/figure_4.png differ
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RMS_false_ASPECT_pre_30_changes.png b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RMS_false_ASPECT_pre_30_changes.png
new file mode 100644
index 00000000000..66515614275
Binary files /dev/null and b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RMS_false_ASPECT_pre_30_changes.png differ
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_25.png b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_25.png
new file mode 100644
index 00000000000..7a49be8d823
Binary files /dev/null and b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_25.png differ
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30.png b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30.png
new file mode 100644
index 00000000000..0a4f717740e
Binary files /dev/null and b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30.png differ
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30_AMG.png b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30_AMG.png
new file mode 100644
index 00000000000..7a585a7589c
Binary files /dev/null and b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30_AMG.png differ
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30_AMG_NSP_none.png b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30_AMG_NSP_none.png
new file mode 100644
index 00000000000..ba2c50e2f70
Binary files /dev/null and b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30_AMG_NSP_none.png differ
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30_AMG_NSP_none_SF_099.png b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30_AMG_NSP_none_SF_099.png
new file mode 100644
index 00000000000..40cc74ed576
Binary files /dev/null and b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_false_ASPECT_30_AMG_NSP_none_SF_099.png differ
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_true_ASPECT_30.png b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_true_ASPECT_30.png
new file mode 100644
index 00000000000..8db1cad9865
Binary files /dev/null and b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_RSM_true_ASPECT_30.png differ
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_iterated_advection_and_stokes.png b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_iterated_advection_and_stokes.png
new file mode 100644
index 00000000000..0a4f717740e
Binary files /dev/null and b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/figure_4_iterated_advection_and_stokes.png differ
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/image_collections/.gitignore b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/image_collections/.gitignore
deleted file mode 100644
index 5e7d2734cfc..00000000000
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/image_collections/.gitignore
+++ /dev/null
@@ -1,4 +0,0 @@
-# Ignore everything in this directory
-*
-# Except this file
-!.gitignore
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/image_data/.gitignore b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/image_data/.gitignore
deleted file mode 100644
index 5e7d2734cfc..00000000000
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/image_data/.gitignore
+++ /dev/null
@@ -1,4 +0,0 @@
-# Ignore everything in this directory
-*
-# Except this file
-!.gitignore
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/image_plotfiles/.gitignore b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/image_plotfiles/.gitignore
deleted file mode 100644
index 5e7d2734cfc..00000000000
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/image_plotfiles/.gitignore
+++ /dev/null
@@ -1,4 +0,0 @@
-# Ignore everything in this directory
-*
-# Except this file
-!.gitignore
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/input.prm b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/input.prm
index a35e871b3ca..7db3e5d4fcb 100644
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/input.prm
+++ b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/input.prm
@@ -1,21 +1,25 @@
 set Additional shared libraries            = ./libdrucker_prager_compositions.so
 set Use years in output instead of seconds = false
-set Nonlinear solver scheme                = iterated Advection and Newton Stokes
+set Nonlinear solver scheme                = single Advection, iterated Newton Stokes
 set Max nonlinear iterations               = 150
 set Nonlinear solver tolerance             = 1e-14
 set Pressure normalization                 = no
 
+# This is an instantaneous benchmark
+set End time                               = 0.0
+
+
 subsection Solver parameters
   subsection Stokes solver parameters
-    set Number of cheap Stokes solver steps = 0
-    set Linear solver tolerance             = 1e-5
+    set Linear solver tolerance = 1e-8
+    set Stokes solver type =  block AMG
   end
 
   subsection Newton solver parameters
     set Max pre-Newton nonlinear iterations      = 5
     set SPD safety factor                        = 0.9
-    set Nonlinear Newton solver switch tolerance = 1e-20 #
-    set Max Newton line search iterations        = 0
+    set Nonlinear Newton solver switch tolerance = 1e-20
+    set Max Newton line search iterations        = 5
     set Maximum linear Stokes solver tolerance   = 9e-1
     set Use Newton residual scaling method       = false
     set Use Newton failsafe                      = false
@@ -94,14 +98,14 @@ subsection Material model
     set Viscosity averaging p       = -1
     set Use deviator of strain-rate = true
     set Use analytical derivative   = false
-    set Reference viscosity         = 5e24
+    set Reference viscosity         = 1e23
     set Reference compressibility   = 0
   end
 end
 
 subsection Mesh refinement
   set Initial adaptive refinement = 0
-  set Initial global refinement   = 5
+  set Initial global refinement   = 4
   set Strategy                    = strain rate
 end
 
@@ -109,6 +113,15 @@ subsection Boundary temperature model
   set Fixed temperature boundary indicators   =
 end
 
+subsection Boundary composition model
+  set Fixed composition boundary indicators = left, right, top
+  set List of model names = function
+
+  subsection Function
+    set Function expression = 0.0
+  end
+end
+
 subsection Boundary velocity model
   # Prescribe a horizontal traction on the vertical boundaries
   set Tangential velocity boundary indicators = 2
@@ -128,7 +141,7 @@ subsection Boundary traction model
 end
 
 subsection Postprocess
-  set List of postprocessors = velocity statistics, pressure statistics, mass flux statistics, visualization
+  set List of postprocessors = velocity statistics, pressure statistics, mass flux statistics, composition statistics, visualization
 
   subsection Visualization
     set List of output variables = material properties,strain rate, spd factor
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/metabash.sh b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/metabash.sh
deleted file mode 100755
index 4932e503272..00000000000
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/metabash.sh
+++ /dev/null
@@ -1,75 +0,0 @@
-#!/bin/bash
-version="c3"
-declare -a grid=(4)
-declare -a NSP=("none") # "symmetric" "none") #Use Newton stabilization preconditioner
-declare -a NSA=("none") # "symmetric" "none") #Use Newton stabilization A block
-declare -a UFS=("truefalse") #"false") # use failsafe
-declare -a RSM=("truefalse") # "true") # use residual scaling
-declare -a LS=(0-5-10) #(10) # line search
-declare -a COMP=(0) # copmressibiltiy
-declare -a phi=(30)
-declare -a ST=("1e-20") # "1e-4")
-declare -a UDS=("false" "true") 
-declare -a SF=("9e-1") # "9.9e-1") # "9.99e-1" "1")
-declare -a AEW=("false") #"true" #"false"
-processes=4
-
-for i_grid in "${grid[@]}"
-do
- for i_NSP in "${NSP[@]}"
- do
-  for i_NSA in "${NSA[@]}"
-  do
-   for i_LS in "${LS[@]}"
-   do
-    for i_RSM in "${RSM[@]}"
-    do
-     for i_UFS in "${UFS[@]}"
-     do
-      for i_COMP in "${COMP[@]}"
-      do
-       for i_phi in "${phi[@]}"
-       do
-        for i_ST in "${ST[@]}"
-        do
-         for i_UDS in "${UDS[@]}"
-         do
-          for i_SF in "${SF[@]}"
-          do
-           for i_AEW in "${AEW[@]}"
-           do
-           output_name="$version""_g""$i_grid""_UFS""$i_UFS""_NSP-""$i_NSP""_NSA-""$i_NSA""_LS""$i_LS""_RSM""$i_RSM""_phi""$i_phi""_COMP""$i_COMP""_ST""$i_ST""_UDS""$i_UDS""_SF""$i_SF""_AEW""$i_AEW"
-           log_file="auto_logs/""$output_name"".log"
-           bash_auto_file="bash_autos/bash_""$output_name"".sh"
-           echo "$output_name"
-           sed \
-            -e "s/declare -a grid.*/declare -a grid=($i_grid)/g" \
-            -e "s/declare -a UFS.*/declare -a UFS=(\"false\" \"true\")/g" \
-            -e "s/declare -a NSP.*/declare -a NSP=($i_NSP)/g" \
-            -e "s/declare -a NSA.*/declare -a NSA=($i_NSA)/g" \
-            -e "s/declare -a ST.*/declare -a ST=($i_ST)/g" \
-            -e "s/declare -a UDS.*/declare -a UDS=($i_UDS)/g" \
-            -e "s/declare -a SF.*/declare -a SF=($i_SF)/g" \
-            -e "s/declare -a AEW.*/declare -a AEW=($i_AEW)/g" \
-            -e "s/declare -a LS.*/declare -a LS=(0 5 10)/g" \
-            -e "s/version=.*/declare -a version=\"$version\"/g" \
-            -e "s/declare -a RSM.*/declare -a RSM=(\"false\" \"true\")/g" \
-            -e "s/phi=.*/phi=\"$i_phi\"/g" \
-            -e "s/COMP=.*/COMP=\"$i_COMP\"/g" \
-            -e "s/processes=.*/processes=\"$processes\"/g" \
-            bash.sh > $bash_auto_file
-           nohup bash ./$bash_auto_file > $log_file 2>&1 &
-           sleep 2 
-           
-           done
-          done
-         done
-        done
-       done
-      done
-     done
-    done
-   done
-  done
- done
-done
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/metabash_itimpes.sh b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/metabash_itimpes.sh
deleted file mode 100755
index 32893cfe197..00000000000
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/metabash_itimpes.sh
+++ /dev/null
@@ -1,76 +0,0 @@
-#!/bin/bash
-version="c3"
-declare -a grid=(4)
-declare -a NSP=("none") # "symmetric" "none") #Use Newton stabilization preconditioner
-declare -a NSA=("none") # "symmetric" "none") #Use Newton stabilization A block
-declare -a UFS=("truefalse") #"false") # use failsafe
-declare -a RSM=("truefalse") # "true") # use residual scaling
-declare -a LS=(0-5-10) #(10) # line search
-declare -a COMP=(0) # copmressibiltiy
-declare -a phi=(30)
-declare -a ST=("1e-20") # "1e-4")
-declare -a UDS=("false" "true") #"true" #"false"
-declare -a SF=("9e-1") # "9.9e-1") # "9.99e-1" "1")
-declare -a AEW=("false") #"true" #"false"
-processes=4
-
-for i_grid in "${grid[@]}"
-do
- for i_NSP in "${NSP[@]}"
- do
-  for i_NSA in "${NSA[@]}"
-  do
-   for i_LS in "${LS[@]}"
-   do
-    for i_RSM in "${RSM[@]}"
-    do
-     for i_UFS in "${UFS[@]}"
-     do
-      for i_COMP in "${COMP[@]}"
-      do
-       for i_phi in "${phi[@]}"
-       do
-        for i_ST in "${ST[@]}"
-        do
-         for i_UDS in "${UDS[@]}"
-         do
-          for i_SF in "${SF[@]}"
-          do
-           for i_AEW in "${AEW[@]}"
-           do
-           output_name="$version""_g""$i_grid""_UFS""$i_UFS""_NSP-""$i_NSP""_NSA-""$i_NSA""_LS""$i_LS""_RSM""$i_RSM""_phi""$i_phi""_COMP""$i_COMP""_ST""$i_ST""_UDS""$i_UDS""_SF""$i_SF""_AEW""$i_AEW"
-           log_file="auto_logs/""$output_name"".log"
-           bash_auto_file="bash_autos/bash_""$output_name"".sh"
-           echo "$output_name"
-           sed \
-            -e "s/declare -a grid.*/declare -a grid=($i_grid)/g" \
-            -e "s/declare -a UFS.*/declare -a UFS=(\"false\" \"true\")/g" \
-            -e "s/declare -a NSP.*/declare -a NSP=($i_NSP)/g" \
-            -e "s/declare -a NSA.*/declare -a NSA=($i_NSA)/g" \
-            -e "s/declare -a ST.*/declare -a ST=($i_ST)/g" \
-            -e "s/declare -a UDS.*/declare -a UDS=($i_UDS)/g" \
-            -e "s/declare -a SF.*/declare -a SF=($i_SF)/g" \
-            -e "s/declare -a AEW.*/declare -a AEW=($i_AEW)/g" \
-            -e "s/declare -a LS.*/declare -a LS=(0 5 10)/g" \
-            -e "s/version=.*/declare -a version=\"$version\"/g" \
-            -e "s/declare -a RSM.*/declare -a RSM=(\"false\" \"true\")/g" \
-            -e "s/phi=.*/phi=\"$i_phi\"/g" \
-            -e "s/COMP=.*/COMP=\"$i_COMP\"/g" \
-            -e "s/processes=.*/processes=\"$processes\"/g" \
-            bash.sh > $bash_auto_file
-           
-           nohup bash ./$bash_auto_file > $log_file 2>&1 &
-           sleep 2 
-           
-           done
-          done
-         done
-        done
-       done
-      done
-     done
-    done
-   done
-  done
- done
-done
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/plot.gnuplot b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/plot.gnuplot
new file mode 100644
index 00000000000..8c39043e08c
--- /dev/null
+++ b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/plot.gnuplot
@@ -0,0 +1,94 @@
+set terminal pngcairo size 1600,1050 enhanced font 'Verdana,14'
+set output 'figure_4.png'
+set yrange [1e-14:1]
+set xrange [0:50]
+set format y '%g'
+set logscale y
+set multiplot title "Spiegelman benchmark convergence result for refinement 4 (64 by 16 cells)" font 'Verdana, 20'
+set size 0.33,0.49
+set origin 0,0.48
+unset key
+set title "Vel: 2.5cm/yr, \\eta_{ref}: 1e23, mLT: 9e-1"
+plot \
+'results/singleAdvectionanditeratedStokes_vel_25_BV_1e23/plot.dat' u 9:10 w l lw 4 lc rgb 'purple' t 'Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P150_vel_25_BV_1e23/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'purple' t 'DC Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P0_vel_25_BV_1e23/plot.dat' u 9:10 w l lw 4 lc rgb 'black' t '0 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P0_vel_25_BV_1e23/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'black' t '0 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P5_vel_25_BV_1e23/plot.dat' u 9:10 w l lw 4 lc rgb 'blue' t '5 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P5_vel_25_BV_1e23/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'blue' t '5 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P15_vel_25_BV_1e23/plot.dat' u 9:10 w l lw 4 lc rgb 'orange' t '15 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P15_vel_25_BV_1e23/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'orange' t '15 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P25_vel_25_BV_1e23/plot.dat' u 9:10 w l lw 4 lc rgb 'red' t '25 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P25_vel_25_BV_1e23/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'red' t '25 Picard stabilized'
+
+set origin 0.33,0.48
+set title "Vel: 5cm/yr, \\eta_{ref}: 1e24, mLT: 9e-1"
+plot \
+'results/singleAdvectionanditeratedStokes_vel_50_BV_1e24/plot.dat' u 9:10 w l lw 4 lc rgb 'purple' t 'Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P150_vel_50_BV_1e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'purple' t 'DC Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P0_vel_50_BV_1e24/plot.dat' u 9:10 w l lw 4 lc rgb 'black' t '0 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P0_vel_50_BV_1e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'black' t '0 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P5_vel_50_BV_1e24/plot.dat' u 9:10 w l lw 4 lc rgb 'blue' t '5 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P5_vel_50_BV_1e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'blue' t '5 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P15_vel_50_BV_1e24/plot.dat' u 9:10 w l lw 4 lc rgb 'orange' t '15 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P15_vel_50_BV_1e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'orange' t '15 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P25_vel_50_BV_1e24/plot.dat' u 9:10 w l lw 4 lc rgb 'red' t '25 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P25_vel_50_BV_1e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'red' t '25 Picard stabilized'
+
+set origin 0.66,0.48
+set title "Vel: 12.5cm/yr, \\eta_{ref}: 5e24, mLT: 9e-1"
+plot \
+'results/singleAdvectionanditeratedStokes_vel_125_BV_5e24/plot.dat' u 9:10 w l lw 4 lc rgb 'purple' t 'Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P150_vel_125_BV_5e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'purple' t 'DC Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P0_vel_125_BV_5e24/plot.dat' u 9:10 w l lw 4 lc rgb 'black' t '0 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P0_vel_125_BV_5e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'black' t '0 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P5_vel_125_BV_5e24/plot.dat' u 9:10 w l lw 4 lc rgb 'blue' t '5 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P5_vel_125_BV_5e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'blue' t '5 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P15_vel_125_BV_5e24/plot.dat' u 9:10 w l lw 4 lc rgb 'orange' t '15 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P15_vel_125_BV_5e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'orange' t '15 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_9e-1_P25_vel_125_BV_5e24/plot.dat' u 9:10 w l lw 4 lc rgb 'red' t '25 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_9e-1_P25_vel_125_BV_5e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'red' t '25 Picard stabilized'
+
+set origin 0.0,0.0
+set title "Vel: 2.5cm/yr, \\eta_{ref}: 1e23, mLT: 1e-8"
+plot \
+'results/singleAdvectionanditeratedStokes_vel_25_BV_1e23/plot.dat' u 9:10 w l lw 4 lc rgb 'purple' t 'Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P150_vel_25_BV_1e23/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'purple' t 'DC Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P0_vel_25_BV_1e23/plot.dat' u 9:10 w l lw 4 lc rgb 'black' t '0 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P0_vel_25_BV_1e23/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'black' t '0 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P5_vel_25_BV_1e23/plot.dat' u 9:10 w l lw 4 lc rgb 'blue' t '5 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P5_vel_25_BV_1e23/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'blue' t '5 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P15_vel_25_BV_1e23/plot.dat' u 9:10 w l lw 4 lc rgb 'orange' t '15 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P15_vel_25_BV_1e23/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'orange' t '15 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P25_vel_25_BV_1e23/plot.dat' u 9:10 w l lw 4 lc rgb 'red' t '25 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P25_vel_25_BV_1e23/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'red' t '25 Picard stabilized'
+
+set origin 0.33,0.0
+set title "Vel: 5cm/yr, \\eta_{ref}: 1e24, mLT: 1e-8"
+plot \
+'results/singleAdvectionanditeratedStokes_vel_50_BV_1e24/plot.dat' u 9:10 w l lw 4 lc rgb 'purple' t 'Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P150_vel_50_BV_1e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'purple' t 'DC Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P0_vel_50_BV_1e24/plot.dat' u 9:10 w l lw 4 lc rgb 'black' t '0 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P0_vel_50_BV_1e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'black' t '0 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P5_vel_50_BV_1e24/plot.dat' u 9:10 w l lw 4 lc rgb 'blue' t '5 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P5_vel_50_BV_1e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'blue' t '5 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P15_vel_50_BV_1e24/plot.dat' u 9:10 w l lw 4 lc rgb 'orange' t '15 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P15_vel_50_BV_1e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'orange' t '15 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P25_vel_50_BV_1e24/plot.dat' u 9:10 w l lw 4 lc rgb 'red' t '25 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P25_vel_50_BV_1e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'red' t '25 Picard stabilized'
+
+set key bottom right
+set key spacing 0.8
+set origin 0.66,0.0
+set title "Vel: 12.5cm/yr, \\eta_{ref}: 5e24, mLT: 1e-8"
+plot \
+'results/singleAdvectionanditeratedStokes_vel_125_BV_5e24/plot.dat' u 9:10 w l lw 4 lc rgb 'purple' t 'Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P150_vel_125_BV_5e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'purple' t 'DC Picard', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P0_vel_125_BV_5e24/plot.dat' u 9:10 w l lw 4 lc rgb 'black' t '0 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P0_vel_125_BV_5e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'black' t '0 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P5_vel_125_BV_5e24/plot.dat' u 9:10 w l lw 4 lc rgb 'blue' t '5 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P5_vel_125_BV_5e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'blue' t '5 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P15_vel_125_BV_5e24/plot.dat' u 9:10 w l lw 4 lc rgb 'orange' t '15 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P15_vel_125_BV_5e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'orange' t '15 Picard stabilized', \
+'results/singleAdvectionandNewtonStokes_NS_none_mLT_1e-8_P25_vel_125_BV_5e24/plot.dat' u 9:10 w l lw 4 lc rgb 'red' t '25 Picard unstabilized', \
+'results/singleAdvectionandNewtonStokes_NS_SPD_mLT_1e-8_P25_vel_125_BV_5e24/plot.dat' u 9:10 w p pt 1 lw 4 lc rgb 'red' t '25 Picard stabilized'
diff --git a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/run.sh b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/run.sh
new file mode 100755
index 00000000000..f5ee0199fcf
--- /dev/null
+++ b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/run.sh
@@ -0,0 +1,66 @@
+#!/bin/bash
+processes=4
+infile="input.prm"
+outfile="output.log"
+errorfile="error.log"
+plotfile="plot.dat"
+
+for case in 1 2 3; do
+  # velocity in m/s
+  if [ $case == 1 ]; then
+    velocity=25
+    U="7.92219116e-11"
+    background_viscosity="1e23"
+  elif [ $case == 2 ]; then
+    velocity=50
+    U="1.58443823e-10"
+    background_viscosity="1e24"
+  elif [ $case == 3 ]; then
+    velocity=125
+    U="3.96109558e-10"
+    background_viscosity="5e24"
+  fi
+
+  # First run one model with only Picard iterations
+  dirname_base="singleAdvectionanditeratedStokes_vel_${velocity}_BV_${background_viscosity}"
+  dirname="results/${dirname_base}"
+  echo "Starting $dirname"
+  mkdir -p $dirname
+
+  sed  \
+    -e "s/set Function expression = if(x<60e3,.*/    set Function expression = if(x<60e3,$U,-$U);0/g" \
+    -e "s/set Reference viscosity .*/    set Reference viscosity = $background_viscosity/g" \
+    -e "s/set Output directory .*/set Output directory = results\/$dirname_base/g" \
+    -e "s/set Nonlinear solver scheme .*/set Nonlinear solver scheme = single Advection, iterated Stokes/g" \
+    -e "s/set List of output variables .*/set List of output variables = material properties,strain rate/g" \
+    -e "s/set Linear solver tolerance .*/set Linear solver tolerance = 1e-14/g" \
+    $infile > "$dirname/$infile"
+
+  nohup mpirun -np $processes ./aspect "${dirname}/${infile}" > ${dirname}/${outfile} 2>${dirname}/${errorfile}
+  grep "Relative nonlinear residual " ${dirname}/${outfile} > ${dirname}/${plotfile}
+
+  # Now run Newton solver models with varying solver parameters
+  for stabilization in "none" "SPD"; do
+    for picard_iterations in 0 5 15 25 150; do # note 150 represents pure defect-correction Picard
+      for maximum_linear_tolerance in "9e-1" "1e-8"; do
+        dirname_base="singleAdvectionandNewtonStokes""_NS_${stabilization}_mLT_${maximum_linear_tolerance}_P${picard_iterations}_vel_${velocity}_BV_${background_viscosity}"
+        dirname="results/$dirname_base"
+        echo "Starting $dirname"
+        mkdir -p $dirname
+
+        sed  \
+          -e "s/set Function expression = if(x<60e3,.*/    set Function expression = if(x<60e3,$U,-$U);0/g" \
+          -e "s/set Reference viscosity .*/    set Reference viscosity = $background_viscosity/g" \
+          -e "s/set Output directory .*/set Output directory = results\/$dirname_base/g" \
+          -e "s/set Max pre-Newton nonlinear iterations .*/    set Max pre-Newton nonlinear iterations = $picard_iterations/g" \
+          -e "s/set Stabilization preconditioner .*/set Stabilization preconditioner = $stabilization/g" \
+          -e "s/set Stabilization velocity block .*/set Stabilization velocity block = $stabilization/g" \
+          -e "s/set Maximum linear Stokes solver tolerance .*/set Maximum linear Stokes solver tolerance = $maximum_linear_tolerance/g" \
+          $infile > "$dirname/$infile"
+
+        nohup mpirun -np $processes ./aspect $dirname/$infile > $dirname/$outfile 2>$dirname/$errorfile
+        grep "Relative nonlinear residual " $dirname/$outfile > $dirname/$plotfile
+      done
+    done
+  done
+done
diff --git a/doc/sphinx/references.bib b/doc/sphinx/references.bib
index c1bf6582649..d7e873fdeb1 100644
--- a/doc/sphinx/references.bib
+++ b/doc/sphinx/references.bib
@@ -12306,4 +12306,15 @@ @article{mcdonough:sun:1995
   pages={223--253},
   year={1995},
   publisher={Elsevier}
-}
\ No newline at end of file
+}
+
+@article{spiegelman:etal:2016,
+  title={On the solvability of incompressible Stokes with viscoplastic rheologies in geodynamics},
+  author={Spiegelman, Marc and May, Dave A and Wilson, Cian R},
+  journal={Geochemistry, Geophysics, Geosystems},
+  volume={17},
+  number={6},
+  pages={2213--2238},
+  year={2016},
+  publisher={Wiley Online Library}
+}