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 b82929c0935..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
@@ -13,7 +13,7 @@ After you have executed `run.sh`, you can use the gnuplot plotting script
 ```{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.
+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/run.sh b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/run.sh
index f5ee0199fcf..a6eee8eed5a 100755
--- a/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/run.sh
+++ b/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/run.sh
@@ -28,8 +28,8 @@ for case in 1 2 3; do
   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 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" \
@@ -49,10 +49,10 @@ for case in 1 2 3; do
         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 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 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" \
diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc
index 6528a67cd56..797d87642d9 100644
--- a/source/simulator/assemblers/newton_stokes.cc
+++ b/source/simulator/assemblers/newton_stokes.cc
@@ -178,6 +178,10 @@ namespace aspect
               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];
@@ -479,6 +483,10 @@ 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];
diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc
index 89c2af8f207..26a034f3ea7 100644
--- a/source/simulator/stokes_matrix_free.cc
+++ b/source/simulator/stokes_matrix_free.cc
@@ -1424,6 +1424,10 @@ namespace aspect
 
                   for (unsigned int q=0; q<n_q_points; ++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];