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..f33a60d10a0 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 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/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..a6eee8eed5a
--- /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}
+}
diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc
index 720260d35ea..797d87642d9 100644
--- a/source/simulator/assemblers/newton_stokes.cc
+++ b/source/simulator/assemblers/newton_stokes.cc
@@ -175,12 +175,19 @@ namespace aspect
           else
             {
               const SymmetricTensor<2,dim> viscosity_derivative_wrt_strain_rate = derivatives->viscosity_derivative_wrt_strain_rate[q];
-              const SymmetricTensor<2,dim> effective_strain_rate =
-                elastic_out == nullptr ? deviator(scratch.material_model_inputs.strain_rate[q]) : elastic_out->viscoelastic_strain_rate[q];
-
               const typename Newton::Parameters::Stabilization
               preconditioner_stabilization = this->get_newton_handler().parameters.preconditioner_stabilization;
 
+              // use the correct strain rate for the Jacobian
+              // when elasticity is enabled use viscoelastic strain rate
+              // when stabilization is enabled, use the deviatoric strain rate because the SPD factor
+              // that is computed is only safe for the deviatoric strain rate (see PR #5580 and issue #5555)
+              SymmetricTensor<2,dim> effective_strain_rate = scratch.material_model_inputs.strain_rate[q];
+              if (elastic_out != nullptr)
+                effective_strain_rate = elastic_out->viscoelastic_strain_rate[q];
+              else if ((preconditioner_stabilization & Newton::Parameters::Stabilization::PD) != Newton::Parameters::Stabilization::none)
+                effective_strain_rate = deviator(effective_strain_rate);
+
               // use the spd factor when the stabilization is PD or SPD
               const double alpha = (preconditioner_stabilization & Newton::Parameters::Stabilization::PD) != Newton::Parameters::Stabilization::none ?
                                    Utilities::compute_spd_factor<dim>(eta, effective_strain_rate, viscosity_derivative_wrt_strain_rate,
@@ -400,9 +407,6 @@ namespace aspect
           const double pressure = scratch.material_model_inputs.pressure[q];
           const double velocity_divergence = scratch.velocity_divergence[q];
 
-          const SymmetricTensor<2,dim> effective_strain_rate =
-            elastic_out == nullptr ? deviator(scratch.material_model_inputs.strain_rate[q]) : elastic_out->viscoelastic_strain_rate[q];
-
           const Tensor<1,dim>
           gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q));
 
@@ -479,6 +483,16 @@ namespace aspect
                   const Newton::Parameters::Stabilization velocity_block_stabilization
                     = this->get_newton_handler().parameters.velocity_block_stabilization;
 
+                  // use the correct strain rate for the Jacobian
+                  // when elasticity is enabled use viscoelastic strain rate
+                  // when stabilization is enabled, use the deviatoric strain rate because the SPD factor
+                  // that is computed is only safe for the deviatoric strain rate (see PR #5580 and issue #5555)
+                  SymmetricTensor<2,dim> effective_strain_rate = scratch.material_model_inputs.strain_rate[q];
+                  if (elastic_out != nullptr)
+                    effective_strain_rate = elastic_out->viscoelastic_strain_rate[q];
+                  else if ((velocity_block_stabilization & Newton::Parameters::Stabilization::PD) != Newton::Parameters::Stabilization::none)
+                    effective_strain_rate = deviator(effective_strain_rate);
+
                   // use the spd factor when the stabilization is PD or SPD
                   const double alpha =  (velocity_block_stabilization & Newton::Parameters::Stabilization::PD)
                                         != Newton::Parameters::Stabilization::none
diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc
index 1a1f3aac534..26a034f3ea7 100644
--- a/source/simulator/stokes_matrix_free.cc
+++ b/source/simulator/stokes_matrix_free.cc
@@ -1424,8 +1424,15 @@ namespace aspect
 
                   for (unsigned int q=0; q<n_q_points; ++q)
                     {
-                      const SymmetricTensor<2,dim> effective_strain_rate =
-                        elastic_out == nullptr ? deviator(in.strain_rate[q]) : elastic_out->viscoelastic_strain_rate[q];
+                      // use the correct strain rate for the Jacobian
+                      // when elasticity is enabled use viscoelastic strain rate
+                      // when stabilization is enabled, use the deviatoric strain rate because the SPD factor
+                      // that is computed is only safe for the deviatoric strain rate (see PR #5580 and issue #5555)
+                      SymmetricTensor<2,dim> effective_strain_rate = in.strain_rate[q];
+                      if (elastic_out != nullptr)
+                        effective_strain_rate = elastic_out->viscoelastic_strain_rate[q];
+                      else if ((sim.newton_handler->parameters.velocity_block_stabilization & Newton::Parameters::Stabilization::PD) != Newton::Parameters::Stabilization::none)
+                        effective_strain_rate = deviator(effective_strain_rate);
 
                       // use the spd factor when the stabilization is PD or SPD.
                       const double alpha =  (sim.newton_handler->parameters.velocity_block_stabilization
diff --git a/tests/spiegelman_et_al_2016.cc b/tests/spiegelman_et_al_2016.cc
new file mode 100644
index 00000000000..7c66b22d012
--- /dev/null
+++ b/tests/spiegelman_et_al_2016.cc
@@ -0,0 +1,21 @@
+/*
+  Copyright (C) 2024 - by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+#include "../benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/drucker_prager_compositions.cc"
diff --git a/tests/spiegelman_et_al_2016.prm b/tests/spiegelman_et_al_2016.prm
new file mode 100644
index 00000000000..5ba18cecdce
--- /dev/null
+++ b/tests/spiegelman_et_al_2016.prm
@@ -0,0 +1,8 @@
+# This is a simple case of the Spiegelman et al. 2016 benchmark set.
+# This test exists to make sure the benchmark input file is executable and
+# the results remain approximately unchanged. See the benchmark for
+# more details on the model and results.
+
+set Dimension = 2
+
+include $ASPECT_SOURCE_DIR/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/input.prm
diff --git a/tests/spiegelman_et_al_2016/screen-output b/tests/spiegelman_et_al_2016/screen-output
new file mode 100644
index 00000000000..ca40c4d0a1f
--- /dev/null
+++ b/tests/spiegelman_et_al_2016/screen-output
@@ -0,0 +1,115 @@
+
+Loading shared library <./libspiegelman_et_al_2016.debug.so>
+
+Number of active cells: 1,024 (on 5 levels)
+Number of degrees of freedom: 18,133 (8,514+1,105+4,257+4,257)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Skipping temperature solve because RHS is zero.
+   Solving C_1 system ... 0 iterations.
+   Initial Newton Stokes residual = 1.56e+14, v = 1.54916e+14, p = 1.83559e+13
+
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 73+0 iterations.
+      Newton system information: Norm of the rhs: 1.56e+14
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 1
+
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 73+0 iterations.
+      Newton system information: Norm of the rhs: 2.96687e+11, Derivative scaling factor: 0
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 2: 0.00190184
+
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 75+0 iterations.
+      Newton system information: Norm of the rhs: 6.45169e+10, Derivative scaling factor: 0
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 3: 0.00041357
+
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 75+0 iterations.
+      Newton system information: Norm of the rhs: 1.68293e+10, Derivative scaling factor: 0
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 4: 0.00010788
+
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 77+0 iterations.
+      Newton system information: Norm of the rhs: 5.77031e+09, Derivative scaling factor: 0
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 5: 3.69892e-05
+
+   Switching from defect correction form of Picard to the Newton solver scheme.
+   The linear solver tolerance is set to 0.342872. Stabilization Preconditioner is SPD and A block is SPD.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 6+0 iterations.
+      Newton system information: Norm of the rhs: 1.90397e+09, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 6: 1.22049e-05
+
+   The linear solver tolerance is set to 0.176941. Stabilization Preconditioner is SPD and A block is SPD.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 14+0 iterations.
+      Newton system information: Norm of the rhs: 3.55866e+08, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 7: 2.28119e-06
+
+   The linear solver tolerance is set to 0.0221025. Stabilization Preconditioner is SPD and A block is SPD.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 21+0 iterations.
+      Newton system information: Norm of the rhs: 4.72273e+07, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 8: 3.02739e-07
+
+   The linear solver tolerance is set to 0.11164. Stabilization Preconditioner is SPD and A block is SPD.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 11+0 iterations.
+      Newton system information: Norm of the rhs: 5.84372e+06, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 9: 3.74598e-08
+
+   The linear solver tolerance is set to 0.0223596. Stabilization Preconditioner is SPD and A block is SPD.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 20+0 iterations.
+      Newton system information: Norm of the rhs: 389976, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 10: 2.49985e-09
+
+   The linear solver tolerance is set to 0.048001. Stabilization Preconditioner is SPD and A block is SPD.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 19+0 iterations.
+      Newton system information: Norm of the rhs: 41104.8, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 11: 2.63493e-10
+
+   The linear solver tolerance is set to 0.0584744. Stabilization Preconditioner is SPD and A block is SPD.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 19+0 iterations.
+      Newton system information: Norm of the rhs: 4468.88, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 12: 2.86467e-11
+
+   The linear solver tolerance is set to 0.0502904. Stabilization Preconditioner is SPD and A block is SPD.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 17+0 iterations.
+      Newton system information: Norm of the rhs: 337.8, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 13: 2.16539e-12
+
+   The linear solver tolerance is set to 0.0403302. Stabilization Preconditioner is SPD and A block is SPD.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 19+0 iterations.
+      Newton system information: Norm of the rhs: 24.3601, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 14: 1.56154e-13
+
+   The linear solver tolerance is set to 0.0418012. Stabilization Preconditioner is SPD and A block is SPD.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 21+0 iterations.
+      Newton system information: Norm of the rhs: 1.87968, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 15: 1.20493e-14
+
+   The linear solver tolerance is set to 0.0435935. Stabilization Preconditioner is SPD and A block is SPD.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 21+0 iterations.
+      Newton system information: Norm of the rhs: 0.231369, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 16: 1.48314e-15
+
+
+   Postprocessing:
+     RMS, max velocity:                  5.21e-11 m/s, 8.72e-11 m/s
+     Pressure min/avg/max:               8.098e+07 Pa, 5.014e+08 Pa, 8.58e+08 Pa
+     Mass fluxes through boundary parts: -0.006483 kg/s, -0.006483 kg/s, 0 kg/s, 0.01297 kg/s
+     Compositions min/max/mass:          0/1/8.705e+08
+     Writing graphical output:           output-spiegelman_et_al_2016/solution/solution-00000
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/spiegelman_et_al_2016_unstabilized.cc b/tests/spiegelman_et_al_2016_unstabilized.cc
new file mode 100644
index 00000000000..7c66b22d012
--- /dev/null
+++ b/tests/spiegelman_et_al_2016_unstabilized.cc
@@ -0,0 +1,21 @@
+/*
+  Copyright (C) 2024 - by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+#include "../benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/drucker_prager_compositions.cc"
diff --git a/tests/spiegelman_et_al_2016_unstabilized.prm b/tests/spiegelman_et_al_2016_unstabilized.prm
new file mode 100644
index 00000000000..798f25750e1
--- /dev/null
+++ b/tests/spiegelman_et_al_2016_unstabilized.prm
@@ -0,0 +1,15 @@
+# Like the test spiegelman_et_al_2016, but without SPD stabilization of
+# the Newton matrix or preconditioner. Because this is the simplest case
+# of the benchmark, convergence is still guaranteed and will be
+# faster than in the stabilized default case.
+
+set Dimension = 2
+
+include $ASPECT_SOURCE_DIR/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/input.prm
+
+subsection Solver parameters
+  subsection Newton solver parameters
+    set Stabilization preconditioner             = none
+    set Stabilization velocity block             = none
+  end
+end
diff --git a/tests/spiegelman_et_al_2016_unstabilized/screen-output b/tests/spiegelman_et_al_2016_unstabilized/screen-output
new file mode 100644
index 00000000000..cf85cf7f436
--- /dev/null
+++ b/tests/spiegelman_et_al_2016_unstabilized/screen-output
@@ -0,0 +1,79 @@
+
+Loading shared library <./libspiegelman_et_al_2016_unstabilized.debug.so>
+
+Number of active cells: 1,024 (on 5 levels)
+Number of degrees of freedom: 18,133 (8,514+1,105+4,257+4,257)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Skipping temperature solve because RHS is zero.
+   Solving C_1 system ... 0 iterations.
+   Initial Newton Stokes residual = 1.56e+14, v = 1.54916e+14, p = 1.83559e+13
+
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 73+0 iterations.
+      Newton system information: Norm of the rhs: 1.56e+14
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 1
+
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 73+0 iterations.
+      Newton system information: Norm of the rhs: 2.96687e+11, Derivative scaling factor: 0
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 2: 0.00190184
+
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 75+0 iterations.
+      Newton system information: Norm of the rhs: 6.45169e+10, Derivative scaling factor: 0
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 3: 0.00041357
+
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 75+0 iterations.
+      Newton system information: Norm of the rhs: 1.68293e+10, Derivative scaling factor: 0
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 4: 0.00010788
+
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 77+0 iterations.
+      Newton system information: Norm of the rhs: 5.77031e+09, Derivative scaling factor: 0
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 5: 3.69892e-05
+
+   Switching from defect correction form of Picard to the Newton solver scheme.
+   The linear solver tolerance is set to 0.342872. Stabilization Preconditioner is none and A block is none.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 6+0 iterations.
+      Newton system information: Norm of the rhs: 1.84613e+09, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 6: 1.18342e-05
+
+   The linear solver tolerance is set to 0.176941. Stabilization Preconditioner is none and A block is none.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 14+0 iterations.
+      Newton system information: Norm of the rhs: 3.08962e+08, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 7: 1.98053e-06
+
+   The linear solver tolerance is set to 5.66517e-05. Stabilization Preconditioner is none and A block is none.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 51+0 iterations.
+      Newton system information: Norm of the rhs: 45303.2, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 8: 2.90405e-10
+
+   The linear solver tolerance is set to 9.48821e-05. Stabilization Preconditioner is none and A block is none.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 53+0 iterations.
+      Newton system information: Norm of the rhs: 3.63378, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 9: 2.32935e-14
+
+   The linear solver tolerance is set to 2.12357e-08. Stabilization Preconditioner is none and A block is none.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 93+0 iterations.
+      Newton system information: Norm of the rhs: 0.0775057, Derivative scaling factor: 1
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 10: 4.96832e-16
+
+
+   Postprocessing:
+     RMS, max velocity:                  5.21e-11 m/s, 8.72e-11 m/s
+     Pressure min/avg/max:               8.098e+07 Pa, 5.014e+08 Pa, 8.58e+08 Pa
+     Mass fluxes through boundary parts: -0.006483 kg/s, -0.006483 kg/s, 0 kg/s, 0.01297 kg/s
+     Compositions min/max/mass:          0/1/8.705e+08
+     Writing graphical output:           output-spiegelman_et_al_2016_unstabilized/solution/solution-00000
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/spiegelman_fail_test/screen-output b/tests/spiegelman_fail_test/screen-output
index 2bccf958fbf..b63ab3f3a74 100644
--- a/tests/spiegelman_fail_test/screen-output
+++ b/tests/spiegelman_fail_test/screen-output
@@ -29,75 +29,75 @@ Number of degrees of freedom: 18,133 (8,514+1,105+4,257+4,257)
    The linear solver tolerance is set to 0.843263. Stabilization Preconditioner is none and A block is none.
    Rebuilding Stokes preconditioner...
    Solving Stokes system (AMG)... 0+1 iterations.
-      Newton system information: Norm of the rhs: 9.54798e+14, Derivative scaling factor: 1
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 1.33204e-16, 0.12306
-      Relative nonlinear residual (total system) after nonlinear iteration 3: 0.12306
+      Newton system information: Norm of the rhs: 9.16397e+14, Derivative scaling factor: 1
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 1.33204e-16, 0.118111
+      Relative nonlinear residual (total system) after nonlinear iteration 3: 0.118111
 
    Skipping temperature solve because RHS is zero.
    Solving C_1 system ... 0 iterations.
    The linear solver tolerance is set to 0.758936. Stabilization Preconditioner is none and A block is none.
    Rebuilding Stokes preconditioner...
    Solving Stokes system (AMG)... 0+1 iterations.
-      Newton system information: Norm of the rhs: 4.36974e+14, Derivative scaling factor: 1
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 1.33204e-16, 0.0563199
-      Relative nonlinear residual (total system) after nonlinear iteration 4: 0.0563199
+      Newton system information: Norm of the rhs: 5.68132e+14, Derivative scaling factor: 1
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 1.33204e-16, 0.0732243
+      Relative nonlinear residual (total system) after nonlinear iteration 4: 0.0732243
 
 
 WARNING: The nonlinear solver in the current timestep failed to converge.
 Acting according to the parameter 'Nonlinear solver failure strategy'...
 Continuing to the next timestep even though solution is not fully converged.
    Postprocessing:
-     RMS, max velocity:                  6.76e-09 m/s, 1.5e-08 m/s
-     Pressure min/avg/max:               -7.019e+10 Pa, 1.063e+09 Pa, 2.587e+11 Pa
-     Mass fluxes through boundary parts: -0.03241 kg/s, -0.03241 kg/s, 0 kg/s, -1.728 kg/s
+     RMS, max velocity:                  2.14e-08 m/s, 2.44e-07 m/s
+     Pressure min/avg/max:               -7.28e+10 Pa, 1.519e+09 Pa, 2.899e+11 Pa
+     Mass fluxes through boundary parts: -0.03241 kg/s, -0.03241 kg/s, 0 kg/s, 4.106 kg/s
      Writing graphical output:           output-spiegelman_fail_test/solution/solution-00000
 
 *** Timestep 1:  t=1 seconds, dt=1 seconds
    Skipping temperature solve because RHS is zero.
-   Solving C_1 system ... 0 iterations.
-   Initial Newton Stokes residual = 8.60842e+15, v = 8.60842e+15, p = 3.41939e+12
+   Solving C_1 system ... 1 iterations.
+   Initial Newton Stokes residual = 8.84992e+15, v = 8.84992e+15, p = 2.98258e+12
 
    Rebuilding Stokes preconditioner...
    Solving Stokes system (AMG)... 0+0 iterations.
-      Newton system information: Norm of the rhs: 4.36222e+14
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 7.2409e-13, 0.0506738
-      Relative nonlinear residual (total system) after nonlinear iteration 1: 0.0506738
+      Newton system information: Norm of the rhs: 5.52475e+14
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 2.43828e-12, 0.0624271
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 0.0624271
 
    Skipping temperature solve because RHS is zero.
    Solving C_1 system ... 0 iterations.
    Rebuilding Stokes preconditioner...
    Solving Stokes system (AMG)... 0+1 iterations.
-      Newton system information: Norm of the rhs: 2.17689e+14, Derivative scaling factor: 0
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 7.2409e-13, 0.025288
-      Relative nonlinear residual (total system) after nonlinear iteration 2: 0.025288
+      Newton system information: Norm of the rhs: 2.2942e+14, Derivative scaling factor: 0
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 1.2658e-13, 0.0259234
+      Relative nonlinear residual (total system) after nonlinear iteration 2: 0.0259234
 
    Skipping temperature solve because RHS is zero.
-   Solving C_1 system ... 0 iterations.
+   Solving C_1 system ... 1 iterations.
    Switching from defect correction form of Picard to the Newton solver scheme.
    The linear solver tolerance is set to 0.843263. Stabilization Preconditioner is none and A block is none.
    Rebuilding Stokes preconditioner...
-   Solving Stokes system (AMG)... 0+1 iterations.
-      Newton system information: Norm of the rhs: 1.35994e+14, Derivative scaling factor: 1
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 6.76328e-13, 0.0157978
-      Relative nonlinear residual (total system) after nonlinear iteration 3: 0.0157978
+   Solving Stokes system (AMG)... 0+3 iterations.
+      Newton system information: Norm of the rhs: 1.92576e+14, Derivative scaling factor: 1
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 1.59356e-12, 0.0217602
+      Relative nonlinear residual (total system) after nonlinear iteration 3: 0.0217602
 
    Skipping temperature solve because RHS is zero.
    Solving C_1 system ... 0 iterations.
    The linear solver tolerance is set to 0.758936. Stabilization Preconditioner is none and A block is none.
    Rebuilding Stokes preconditioner...
    Solving Stokes system (AMG)... 0+1 iterations.
-      Newton system information: Norm of the rhs: 7.18526e+13, Derivative scaling factor: 1
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 5.87227e-13, 0.00834678
-      Relative nonlinear residual (total system) after nonlinear iteration 4: 0.00834678
+      Newton system information: Norm of the rhs: 1.08722e+14, Derivative scaling factor: 1
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 2.92488e-13, 0.0122851
+      Relative nonlinear residual (total system) after nonlinear iteration 4: 0.0122851
 
 
 WARNING: The nonlinear solver in the current timestep failed to converge.
 Acting according to the parameter 'Nonlinear solver failure strategy'...
 Continuing to the next timestep even though solution is not fully converged.
    Postprocessing:
-     RMS, max velocity:                  2.71e-09 m/s, 8.03e-09 m/s
-     Pressure min/avg/max:               -2.037e+10 Pa, 1.602e+08 Pa, 3.381e+10 Pa
-     Mass fluxes through boundary parts: -0.03241 kg/s, -0.03241 kg/s, 0 kg/s, -1.176 kg/s
+     RMS, max velocity:                  2.94e-09 m/s, 1.54e-08 m/s
+     Pressure min/avg/max:               -2.771e+10 Pa, 3.617e+08 Pa, 6.167e+10 Pa
+     Mass fluxes through boundary parts: -0.03241 kg/s, -0.03241 kg/s, 0 kg/s, -0.1108 kg/s
 
 Termination requested by criterion: end time