diff --git a/benchmarks/free_surface_tractions/viscoelastic/README b/benchmarks/free_surface_tractions/viscoelastic/README
index f2043b04626..468045f5cf4 100644
--- a/benchmarks/free_surface_tractions/viscoelastic/README
+++ b/benchmarks/free_surface_tractions/viscoelastic/README
@@ -11,7 +11,9 @@ The benchmark compares their solution for the surface displacement
 
 A surface pressure of rho_l*g*H0 (where rho_l is the load density,
  H0 is the load height) is applied on the surface for r<r0 (where
- r0 is the load radius), for t>0. The load is fully removed by t=t1,
+ r0 is the load radius), for t>0. The load is either instantaneously
+ loaded/emplaced at times t0<t<t1, or thinned linearly over this
+ time interval ("linear unloading"). The load is fully removed by t=t1,
  and the surface relaxes. This is done both in a 2-D and 3-D geometry
  (by symmetry the load is centered on the left boundary or left/front
  corner). The input files are:
@@ -22,7 +24,7 @@ A surface pressure of rho_l*g*H0 (where rho_l is the load density,
 
 Changing the stress instantaneously as the load is applied/removed
  leads to oscillations in the elastic response. This can be addressed
- by setting a fixed elastic timestep larger than the numerical 
+ by setting a fixed elastic timestep equal to the numerical ("Maximum")
  timestep and turning on stress averaging:
    'free_surface_VE_cylinder_2D_loading_fixed_elastic_dt.prm'
    'free_surface_VE_cylinder_3D_loading_fixed_elastic_dt.prm'
@@ -36,4 +38,19 @@ The 'topography' output files may be compared against the analytical
  deflection through time, 'compare_viscous_def_profile.gnuplot' for
  deflection  of profile through time). Equivalent scripts for the
  unloading case are also provided.
+ 
+Note that while the analytical and numerical results for the deflection
+ of the surface agree well near the center of the load (left boundary),
+ the solutions do not match as well on the right (free-slip) boundary.
+ The numerical free surface rides up vertically against the right
+ boundary to conserve mass, while the analytical solution assumes an
+ infinite half-space, predicting near-zero displacement far from the
+ load center. This is also true for the viscous benchmark. This problem
+ is less pronounced in 3-D as the extra mass may be distributed over
+ a larger area. An open far (right/back) boundary resolves this problem
+ in 3-D.
 
+The solutions match well in 3-D. In 2-D, the geometries of the loading
+ function are different (Cartesian in ASPECT vs cylindrical analytically).
+ As such, the agreement in 2-D breaks down for small r0 (load width) or
+ if the right boundary is placed further away.
diff --git a/benchmarks/free_surface_tractions/viscoelastic/compare_VE_def.gnuplot b/benchmarks/free_surface_tractions/viscoelastic/compare_VE_def.gnuplot
index 509aaeb93e7..d9f375c1bed 100644
--- a/benchmarks/free_surface_tractions/viscoelastic/compare_VE_def.gnuplot
+++ b/benchmarks/free_surface_tractions/viscoelastic/compare_VE_def.gnuplot
@@ -6,6 +6,6 @@ set title "Depression of surface in response to loading (inst. from 0 to 1000 ye
 set xlabel "Time [years]"
 set ylabel "Maximum depression of surface [m]"
 
-plot 'output_free_surface_VE_cylinder_2D_loading/statistics' using 2:27 title 'Numerical, no stress averaging' with linesp
-replot 'output_free_surface_VE_cylinder_2D_loading_fixed_elastic_dt/statistics' using 2:27 title 'Numerical, with stress averaging' with linesp
+plot 'output_free_surface_VE_cylinder_2D_loading/statistics' using 2:28 title 'Numerical, no stress averaging' with linesp
+replot 'output_free_surface_VE_cylinder_2D_loading_fixed_elastic_dt/statistics' using 2:28 title 'Numerical, with stress averaging' with linesp
 replot 'soln_z0.txt' using 1:2 title 'Analytical (Nakiboglu and Lambeck, 1982)'
diff --git a/benchmarks/free_surface_tractions/viscoelastic/compare_VE_def_unloading.gnuplot b/benchmarks/free_surface_tractions/viscoelastic/compare_VE_def_unloading.gnuplot
index 3d81080fb3d..a7fe9b7ed4f 100644
--- a/benchmarks/free_surface_tractions/viscoelastic/compare_VE_def_unloading.gnuplot
+++ b/benchmarks/free_surface_tractions/viscoelastic/compare_VE_def_unloading.gnuplot
@@ -6,5 +6,5 @@ set title "Depression of surface in response to linear unloading (from 0 to 1000
 set xlabel "Time [years]"
 set ylabel "Maximum depression of surface [m]"
 
-plot 'output_free_surface_VE_cylinder_2D_loading_unloading/statistics' using 2:27 title 'Numerical, no stress averaging' with linesp
+plot 'output_free_surface_VE_cylinder_2D_loading_unloading/statistics' using 2:28 title 'Numerical, no stress averaging' with linesp
 replot 'soln2_z0.txt' using 1:2 title 'Analytical (Nakiboglu and Lambeck, 1982)'
diff --git a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading.prm b/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading.prm
index da33468c7c2..c20ed873055 100644
--- a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading.prm
+++ b/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading.prm
@@ -50,10 +50,13 @@ subsection Formulation
   set Enable elasticity = true
 end
 
-subsection Free surface
-  set Free surface boundary indicators = top
-  set Surface velocity projection      = normal
+subsection Mesh deformation
   set Additional tangential mesh velocity boundary indicators = left,right
+  set Mesh deformation boundary indicators = top: free surface
+
+  subsection Free surface
+    set Surface velocity projection      = normal
+  end
 end
 
 # Velocity boundary conditions
diff --git a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading_fixed_elastic_dt.prm b/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading_fixed_elastic_dt.prm
index adc005c3661..03a48f144d2 100644
--- a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading_fixed_elastic_dt.prm
+++ b/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading_fixed_elastic_dt.prm
@@ -1,6 +1,6 @@
-# Same as free_surface_VE_cylinder_2D_loading.prm, instead
-# setting fixed elastic timestep and allowing stress averaging,
-# to smooth out elastic response to instantaneous stress changes.
+# Same as free_surface_VE_cylinder_2D_loading.prm, but uses
+# a fixed elastic timestep in combination with stress averaging
+# to smooth out the elastic response to instantaneous stress changes.
 
 include free_surface_VE_cylinder_2D_loading.prm
 
@@ -21,4 +21,3 @@ subsection Material model
   end
 
 end
-
diff --git a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_3D_loading.prm b/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_3D_loading.prm
index 39c4612a582..74b3c473051 100644
--- a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_3D_loading.prm
+++ b/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_3D_loading.prm
@@ -51,10 +51,13 @@ subsection Formulation
   set Enable elasticity = true
 end
 
-subsection Free surface
-  set Free surface boundary indicators = top
-  set Surface velocity projection      = normal
-  set Additional tangential mesh velocity boundary indicators = left,right,front,back
+subsection Mesh deformation
+  set Additional tangential mesh velocity boundary indicators = left, right, front, back
+  set Mesh deformation boundary indicators = top: free surface
+
+  subsection Free surface
+    set Surface velocity projection      = normal
+  end
 end
 
 # Velocity boundary conditions
diff --git a/benchmarks/free_surface_tractions/viscous/README b/benchmarks/free_surface_tractions/viscous/README
index 1f6d534da17..190e6ed9b02 100644
--- a/benchmarks/free_surface_tractions/viscous/README
+++ b/benchmarks/free_surface_tractions/viscous/README
@@ -3,16 +3,19 @@ Benchmark of infinite viscous half-space loaded/unloaded by an
  The Motion of a Viscous Fluid Under a Surface Load (1935), Physics
  6, 265; doi:10.1063/1.1745329
 
+This benchmark is identical to the viscoelastic benchmarks, in the
+ limit of negligible elasticity.
+
 The benchmark compares his solution for the surface displacement
  following the application of an instantaneous cylindrical surface
  load (Eq. 3.34) with the deformation of the free surface in ASPECT
  after imposing surface boundary tractions.
 
 A surface pressure of rho_l*g*H0 (where rho_l is the load density,
- H0 is the load height) is applied on the surface for r<r0 (where
- r0 is the load radius), for t>0. This is done both in a 2-D and 
- 3-D geometry (by symmetry the load is centered on the left boundary
- or left/front corner). The input files are: 
+ H0 is the load height) is applied instantaneously on the surface
+ for r<r0 (where r0 is the load radius), for t>0. This is done both
+ in a 2-D and 3-D geometry (by symmetry the load is centered on the
+ left boundary or left/front corner). The input files are: 
    'free_surface_viscous_cylinder_2D_loading.prm'
    'free_surface_viscous_cylinder_3D_loading.prm'
 
@@ -32,4 +35,20 @@ The 'topography' output files may be compared against an analytical
  ASPECT 'topography' output files, using the provided gnuplot script
  ('compare_viscous_def.gnuplot' for maximum surface deflection
  through time, 'compare_viscous_def_profile.gnuplot' for deflection
- of profile through time). 
+ of profile through time).
+ 
+Note that while the analytical and numerical results for the deflection
+ of the surface agree well near the center of the load (left boundary),
+ the solutions do not match as well on the right (free-slip) boundary.
+ The numerical free surface rides up vertically against the right
+ boundary to conserve mass, while the analytical solution assumes an
+ infinite half-space, predicting near-zero displacement far from the
+ load center. This is also true for the visoelastic benchmark. This problem
+ is less pronounced in 3-D as the extra mass may be distributed over
+ a larger area. An open far (right/back) boundary resolves this problem
+ in 3-D.
+
+The solutions match well in 3-D. In 2-D, the geometries of the loading
+ function are different (Cartesian in ASPECT vs cylindrical analytically).
+ As such, the agreement in 2-D breaks down for small r0 (load width) or
+ if the right boundary is placed further away.
diff --git a/benchmarks/free_surface_tractions/viscous/free_surface_viscous_cylinder_2D_loading.prm b/benchmarks/free_surface_tractions/viscous/free_surface_viscous_cylinder_2D_loading.prm
index 526bbcc8cc9..fe7a868f60a 100644
--- a/benchmarks/free_surface_tractions/viscous/free_surface_viscous_cylinder_2D_loading.prm
+++ b/benchmarks/free_surface_tractions/viscous/free_surface_viscous_cylinder_2D_loading.prm
@@ -5,7 +5,7 @@
 # axisymmetric cylindrical load over a viscous half-space.
 
 # This is done by using the viscoelastic model and 
-# setting a high shear modulus.
+# setting a high shear modulus, which produces viscous behavior. 
 
 include ../viscoelastic/free_surface_VE_cylinder_2D_loading.prm
 
@@ -22,4 +22,3 @@ subsection Material model
   end
 
 end
-