diff --git a/benchmarks/particle_integration_scheme/circle.prm b/benchmarks/particle_integration_scheme/circle.prm
new file mode 100644
index 00000000000..f6713c55a8a
--- /dev/null
+++ b/benchmarks/particle_integration_scheme/circle.prm
@@ -0,0 +1,80 @@
+set Dimension = 2
+set Start time = 0
+set End time = 1.0
+set Use years in output instead of seconds = false
+set CFL number = 1.0
+set Output directory = circle_euler_1.0
+set Timing output frequency = 100
+
+subsection Geometry model
+ set Model name = box
+ subsection Box
+ set X extent = 4
+ set Y extent = 4
+ set Box origin X coordinate = -2
+ set Box origin Y coordinate = -2
+ end
+end
+
+
+subsection Prescribed Stokes solution
+ set Model name = function
+ subsection Velocity function
+ set Variable names = x,y,t
+ set Function expression = 2 * pi * (1+0.0*sin(t))*-y; 2 * pi * (1+0.0*sin(t))*x
+ end
+end
+
+set Nonlinear solver scheme = Advection only
+
+subsection Initial temperature model
+ set Model name = function
+end
+
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 0
+ end
+end
+
+
+subsection Material model
+ set Model name = simple
+end
+
+
+subsection Mesh refinement
+ set Initial global refinement = 4
+ set Initial adaptive refinement = 0
+ set Time steps between mesh refinement = 0
+end
+
+############### Parameters describing what to do with the solution
+
+subsection Postprocess
+ set List of postprocessors = visualization, particles
+
+ subsection Visualization
+ set Output format = vtu
+ set Time between graphical output = 1
+ end
+
+ subsection Particles
+ set Number of particles = 50
+ set Time between data output = 0.1
+ set Data output format = ascii
+
+ set Particle generator name = ascii file
+ set Integration scheme = rk2
+
+ subsection Generator
+ subsection Ascii file
+ set Data directory = $ASPECT_SOURCE_DIR/benchmarks/particle_integration_scheme/
+ set Data file name = particle.dat
+ end
+ end
+ end
+end
diff --git a/benchmarks/particle_integration_scheme/create.sh b/benchmarks/particle_integration_scheme/create.sh
new file mode 100644
index 00000000000..30bef13e0dc
--- /dev/null
+++ b/benchmarks/particle_integration_scheme/create.sh
@@ -0,0 +1,31 @@
+#!/bin/bash
+
+for int in euler rk2 rk4; do #euler rk2 rk4
+ rm max_deviation_$int
+
+ for cfl in 1.0 0.5 0.25 0.125 0.0625 0.03125; do
+ # Prepare and run model
+ OUTPUT_DIR=circle_${int}_${cfl}
+ echo set Output directory = $OUTPUT_DIR > output_dir.prm
+ echo set CFL number = $cfl > cfl.prm
+
+ cat circle.prm output_dir.prm cfl.prm ${int}.prm | aspect --
+
+ # format and print results
+ tail -n +8 $OUTPUT_DIR/particles/particles-00000.0000.gnuplot | sort -g -k 3 > sort_temp_0
+ tail -n +8 $OUTPUT_DIR/particles/particles-00001.0000.gnuplot | sort -g -k 3 > sort_temp_1
+ paste sort_temp_0 sort_temp_1 > sort_temp
+ gawk -v cfl=$cfl 'BEGIN{max_deviation=0.0} {deviation=sqrt(($1-$4)*($1-$4)+($2-$5)*($2-$5));if(deviation>max_deviation) {max_deviation=deviation}} END{print cfl, max_deviation}' sort_temp >> max_deviation_$int
+
+
+ done
+done
+
+# make visualization
+echo set Output directory = circle_vis > output_dir.prm
+cat circle.prm vis.prm output_dir.prm | aspect --
+
+# make plot
+#gnuplot plot_convergence.plt
+
+rm output_dir.prm cfl.prm sort_temp_0 sort_temp_1 sort_temp
diff --git a/benchmarks/particle_integration_scheme/doc/convergence.svg b/benchmarks/particle_integration_scheme/doc/convergence.svg
new file mode 100644
index 00000000000..56e1fcc1c82
--- /dev/null
+++ b/benchmarks/particle_integration_scheme/doc/convergence.svg
@@ -0,0 +1,493 @@
+
+
+
diff --git a/benchmarks/particle_integration_scheme/doc/particle_integration_scheme.md b/benchmarks/particle_integration_scheme/doc/particle_integration_scheme.md
new file mode 100644
index 00000000000..474da9c1579
--- /dev/null
+++ b/benchmarks/particle_integration_scheme/doc/particle_integration_scheme.md
@@ -0,0 +1,27 @@
+# Particle integration scheme benchmark
+
+*This section was contributed by Gabriel Johnston and Rene Gassmöller*
+
+This benchmark is designed to test the convergence of various particle integration schemes (Forward-Euler, RK2, RK4) as detailed in {cite}`gassmoller:etal:2018`. The first benchmark involves 100 particles moving in circular motion around a ring in a 2D box. The error is measured as the distance between the final and initial positions, which should ideally be zero. The second benchmark involves particles advected in a flow field defined by {math:numref}`eq:flow`.
+
+```{math}
+:label: eq:flow
+\mathbf{u}_x(t) = 1 + e^t \quad \text{and} \quad \mathbf{u}_y(t) = 1 + e^t.
+```
+The final position analytic position is represented by {math:numref}`eq:fposition` and the distance between the initial and final position is defined by {math:numref}`eq:fidistance`.
+
+```{math}
+:label: eq:fposition
+\mathbf{x}(1) = \mathbf{x}_0 + \int_{0}^{1} \mathbf{u}(t) \, dt
+```
+
+```{math}
+:label: eq:fidistance
+d = \| \mathbf{x}_1 - \mathbf{x}_0 \|_2 = \sqrt{2} e
+```
+
+```{figure-md} fig:particle-integration-scheme
+
+
+Figure shows the convergence order of particle advection schemes in spatially variable flow (a) and temporally changing flow (b). Additionally, the figure (right) depicts the RK4 scheme using the analytical velocity solution for the halfway time instead of interpolating the discrete time steps. Figure adapted from {cite}`gassmoller:etal:2018`
+```
diff --git a/benchmarks/particle_integration_scheme/euler.prm b/benchmarks/particle_integration_scheme/euler.prm
new file mode 100644
index 00000000000..5670d14da70
--- /dev/null
+++ b/benchmarks/particle_integration_scheme/euler.prm
@@ -0,0 +1,5 @@
+subsection Postprocess
+ subsection Particles
+ set Integration scheme = euler
+ end
+end
diff --git a/benchmarks/particle_integration_scheme/particle.dat b/benchmarks/particle_integration_scheme/particle.dat
new file mode 100644
index 00000000000..b89bb3cb3d1
--- /dev/null
+++ b/benchmarks/particle_integration_scheme/particle.dat
@@ -0,0 +1,100 @@
+1.0000 0.0000
+0.9980 0.0634
+0.9920 0.1266
+0.9819 0.1893
+0.9679 0.2511
+0.9501 0.3120
+0.9284 0.3717
+0.9029 0.4298
+0.8738 0.4862
+0.8413 0.5406
+0.8053 0.5929
+0.7660 0.6428
+0.7237 0.6901
+0.6785 0.7346
+0.6306 0.7761
+0.5801 0.8146
+0.5272 0.8497
+0.4723 0.8815
+0.4154 0.9096
+0.3569 0.9341
+0.2969 0.9549
+0.2358 0.9718
+0.1736 0.9848
+0.1108 0.9938
+0.0476 0.9989
+-0.0159 0.9999
+-0.0792 0.9969
+-0.1423 0.9898
+-0.2048 0.9788
+-0.2665 0.9638
+-0.3271 0.9450
+-0.3863 0.9224
+-0.4441 0.8960
+-0.5000 0.8660
+-0.5539 0.8326
+-0.6056 0.7958
+-0.6549 0.7557
+-0.7015 0.7127
+-0.7453 0.6668
+-0.7861 0.6182
+-0.8237 0.5671
+-0.8580 0.5137
+-0.8888 0.4582
+-0.9161 0.4009
+-0.9397 0.3420
+-0.9595 0.2817
+-0.9754 0.2203
+-0.9874 0.1580
+-0.9955 0.0951
+-0.9995 0.0317
+-0.9995 -0.0317
+-0.9955 -0.0951
+-0.9874 -0.1580
+-0.9754 -0.2203
+-0.9595 -0.2817
+-0.9397 -0.3420
+-0.9161 -0.4009
+-0.8888 -0.4582
+-0.8580 -0.5137
+-0.8237 -0.5671
+-0.7861 -0.6182
+-0.7453 -0.6668
+-0.7015 -0.7127
+-0.6549 -0.7557
+-0.6056 -0.7958
+-0.5539 -0.8326
+-0.5000 -0.8660
+-0.4441 -0.8960
+-0.3863 -0.9224
+-0.3271 -0.9450
+-0.2665 -0.9638
+-0.2048 -0.9788
+-0.1423 -0.9898
+-0.0792 -0.9969
+-0.0159 -0.9999
+0.0476 -0.9989
+0.1108 -0.9938
+0.1736 -0.9848
+0.2358 -0.9718
+0.2969 -0.9549
+0.3569 -0.9341
+0.4154 -0.9096
+0.4723 -0.8815
+0.5272 -0.8497
+0.5801 -0.8146
+0.6306 -0.7761
+0.6785 -0.7346
+0.7237 -0.6901
+0.7660 -0.6428
+0.8053 -0.5929
+0.8413 -0.5406
+0.8738 -0.4862
+0.9029 -0.4298
+0.9284 -0.3717
+0.9501 -0.3120
+0.9679 -0.2511
+0.9819 -0.1893
+0.9920 -0.1266
+0.9980 -0.0634
+1.0000 -0.0000
diff --git a/benchmarks/particle_integration_scheme/plot_convergence.plt b/benchmarks/particle_integration_scheme/plot_convergence.plt
new file mode 100644
index 00000000000..fe0f430e2f4
--- /dev/null
+++ b/benchmarks/particle_integration_scheme/plot_convergence.plt
@@ -0,0 +1,63 @@
+set terminal svg size 1500,900 font "Arial,20" enhanced background rgb 'white'
+set output 'convergence.svg'
+
+set multiplot layout 1,2
+
+# Plot A
+set size ratio 1.0
+set origin 0.0,0.0
+set lmargin at screen 0.1
+set rmargin at screen 0.45
+set tmargin at screen 0.9
+set bmargin at screen 0.25
+
+set xrange[2.5e-2:1.2]
+set yrange[1e-14:1.0]
+set xlabel "CFL Number" font "Arial,18"
+set ylabel "L2 errors" font "Arial,18"
+set logscale xy
+
+f(x) = 0.1 * x
+g(x) = 0.005 * x*x
+h(x) = 1e-6 * x*x*x*x
+set key font "Arial,16" outside center bottom samplen 1 maxrows 3 box
+
+set label 1 "A)" at screen 0.25, screen 0.95 font "Arial,22"
+
+plot "spatial/max_deviation_euler" using 1:2 with linespoints linetype 1 lw 3 pointtype 7 ps 1.0 linecolor rgb "#0A5F02" title "Forward Euler", \
+ "spatial/max_deviation_rk2" using 1:2 with linespoints linetype 1 lw 3 pointtype 9 ps 1.0 linecolor rgb "#1c567a" title "RK2", \
+ "spatial/max_deviation_rk4" using 1:2 with linespoints linetype 1 lw 3 pointtype 11 ps 1.0 linecolor rgb "#814292" title "RK4", \
+ f(x) title 'linear' with lines dashtype 2 lw 3 linecolor rgb "gray", \
+ g(x) title 'quadratic' with lines dashtype 3 lw 3 linecolor rgb "gray", \
+ h(x) title 'quartic' with lines dashtype 4 lw 3 linecolor rgb "gray"
+
+# Plot B
+set size ratio 1.0
+set origin 0.5,0.0
+set lmargin at screen 0.55
+set rmargin at screen 0.9
+set tmargin at screen 0.9
+set bmargin at screen 0.25
+
+set xrange[2.5e-2:1.2]
+set yrange[1e-14:1.0]
+set xlabel "CFL Number" font "Arial,18"
+set ylabel "L2 errors" font "Arial,18"
+set logscale xy
+
+f(x) = 0.05 * x
+g(x) = 2e-3 * x*x
+h(x) = 1e-7 * x*x*x*x
+set key font "Arial,16" outside center bottom samplen 1 maxrows 4 box
+
+set label 2 "B)" at screen 0.75, screen 0.95 font "Arial,22"
+
+plot "time_exponential/max_deviation_euler" using 1:2 with linespoints linetype 1 lw 3 pointtype 7 ps 1.0 linecolor rgb "#0A5F02" notitle , \
+ "time_exponential/max_deviation_rk2" using 1:2 with linespoints linetype 1 lw 3 pointtype 8 ps 1.0 linecolor rgb "#1c567a" notitle , \
+ "time_exponential/max_deviation_rk4" using 1:2 with linespoints linetype 1 lw 3 pointtype 11 ps 1.0 linecolor rgb "#814292" notitle, \
+ "time_exponential/max_deviation_rk4_analytic" using 1:2 with linespoints linetype 2 lw 3 pointtype 13 ps 1.0 linecolor rgb "#814292" title "RK4 with analytic halfstep", \
+ f(x) notitle with lines dashtype 2 lw 3 linecolor rgb "gray", \
+ g(x) notitle with lines dashtype 3 lw 3 linecolor rgb "gray", \
+ h(x) notitle with lines dashtype 4 lw 3 linecolor rgb "gray"
+
+unset multiplot
diff --git a/benchmarks/particle_integration_scheme/rk2.prm b/benchmarks/particle_integration_scheme/rk2.prm
new file mode 100644
index 00000000000..8692b18e5e3
--- /dev/null
+++ b/benchmarks/particle_integration_scheme/rk2.prm
@@ -0,0 +1,5 @@
+subsection Postprocess
+ subsection Particles
+ set Integration scheme = rk2
+ end
+end
diff --git a/benchmarks/particle_integration_scheme/rk4.prm b/benchmarks/particle_integration_scheme/rk4.prm
new file mode 100644
index 00000000000..b40939cc7d9
--- /dev/null
+++ b/benchmarks/particle_integration_scheme/rk4.prm
@@ -0,0 +1,5 @@
+subsection Postprocess
+ subsection Particles
+ set Integration scheme = rk4
+ end
+end
diff --git a/doc/modules/changes/20240601_orionjohnston b/doc/modules/changes/20240601_orionjohnston
new file mode 100644
index 00000000000..d96b0846995
--- /dev/null
+++ b/doc/modules/changes/20240601_orionjohnston
@@ -0,0 +1,3 @@
+New: Benchmark from Gassmoller et al. 2018 for particle integration schemes.
+
+(Gabriel Johnston, 2024/06/01)
diff --git a/doc/sphinx/user/benchmarks/index.md b/doc/sphinx/user/benchmarks/index.md
index 1cce19d4b18..ce39b99eec6 100644
--- a/doc/sphinx/user/benchmarks/index.md
+++ b/doc/sphinx/user/benchmarks/index.md
@@ -50,6 +50,7 @@ benchmarks/inclusion/doc/inclusion.md
benchmarks/burstedde/doc/burstedde.md
benchmarks/slab_detachment/doc/slab_detachment.md
benchmarks/hollow_sphere/doc/hollow_sphere.md
+benchmarks/particle_integration_scheme/doc/particle_integration_scheme.md
benchmarks/annulus/doc/annulus.md
benchmarks/rigid_shear/doc/rigid_shear.md
../cookbooks/cookbooks/stokes/doc/stokes.md
diff --git a/tests/benchmark_particle_integration_scheme.prm b/tests/benchmark_particle_integration_scheme.prm
new file mode 100644
index 00000000000..dbff68e5758
--- /dev/null
+++ b/tests/benchmark_particle_integration_scheme.prm
@@ -0,0 +1,7 @@
+# A benchmark of particle integration schemes in Gassmoller 2018
+# See the benchmark folder for more information.
+set Dimension = 2
+
+include $ASPECT_SOURCE_DIR/benchmarks/particle_integration_scheme/circle.prm
+
+set End time = 0.1
diff --git a/tests/benchmark_particle_integration_scheme/particles/particles-00000.0000.gnuplot b/tests/benchmark_particle_integration_scheme/particles/particles-00000.0000.gnuplot
new file mode 100644
index 00000000000..c6579db76c7
--- /dev/null
+++ b/tests/benchmark_particle_integration_scheme/particles/particles-00000.0000.gnuplot
@@ -0,0 +1,207 @@
+# This file was generated by the deal.II library.
+# Date = 2024/6/1
+# Time = 16:23:48
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+#
+-0.6549 -0.7557 63
+
+-0.6056 -0.7958 64
+
+-0.5539 -0.8326 65
+
+-0.5 -0.866 66
+
+-0.858 -0.5137 58
+
+-0.8237 -0.5671 59
+
+-0.7861 -0.6182 60
+
+-0.7453 -0.6668 61
+
+-0.7015 -0.7127 62
+
+-0.4441 -0.896 67
+
+-0.3863 -0.9224 68
+
+-0.3271 -0.945 69
+
+-0.2665 -0.9638 70
+
+-0.2048 -0.9788 71
+
+-0.1423 -0.9898 72
+
+-0.0792 -0.9969 73
+
+-0.0159 -0.9999 74
+
+-0.9595 -0.2817 54
+
+-0.9397 -0.342 55
+
+-0.9161 -0.4009 56
+
+-0.8888 -0.4582 57
+
+-0.9995 -0.0317 50
+
+-0.9955 -0.0951 51
+
+-0.9874 -0.158 52
+
+-0.9754 -0.2203 53
+
+0.0476 -0.9989 75
+
+0.1108 -0.9938 76
+
+0.1736 -0.9848 77
+
+0.2358 -0.9718 78
+
+0.2969 -0.9549 79
+
+0.3569 -0.9341 80
+
+0.4154 -0.9096 81
+
+0.4723 -0.8815 82
+
+0.5272 -0.8497 83
+
+0.5801 -0.8146 84
+
+0.6306 -0.7761 85
+
+0.6785 -0.7346 86
+
+0.7237 -0.6901 87
+
+0.766 -0.6428 88
+
+0.8053 -0.5929 89
+
+0.8413 -0.5406 90
+
+0.8738 -0.4862 91
+
+0.9029 -0.4298 92
+
+0.9284 -0.3717 93
+
+0.9501 -0.312 94
+
+0.9679 -0.2511 95
+
+1 0 0
+
+0.9819 -0.1893 96
+
+0.992 -0.1266 97
+
+0.998 -0.0634 98
+
+1 0 99
+
+-0.9754 0.2203 46
+
+-0.9874 0.158 47
+
+-0.9955 0.0951 48
+
+-0.9995 0.0317 49
+
+-0.8888 0.4582 42
+
+-0.9161 0.4009 43
+
+-0.9397 0.342 44
+
+-0.9595 0.2817 45
+
+-0.7861 0.6182 39
+
+-0.8237 0.5671 40
+
+-0.858 0.5137 41
+
+-0.7015 0.7127 37
+
+-0.7453 0.6668 38
+
+-0.5 0.866 33
+
+-0.5539 0.8326 34
+
+-0.6056 0.7958 35
+
+-0.6549 0.7557 36
+
+-0.2665 0.9638 29
+
+-0.3271 0.945 30
+
+-0.3863 0.9224 31
+
+-0.4441 0.896 32
+
+-0.0159 0.9999 25
+
+-0.0792 0.9969 26
+
+-0.1423 0.9898 27
+
+-0.2048 0.9788 28
+
+0.998 0.0634 1
+
+0.992 0.1266 2
+
+0.9819 0.1893 3
+
+0.9679 0.2511 4
+
+0.9501 0.312 5
+
+0.9284 0.3717 6
+
+0.9029 0.4298 7
+
+0.8738 0.4862 8
+
+0.2358 0.9718 21
+
+0.1736 0.9848 22
+
+0.1108 0.9938 23
+
+0.0476 0.9989 24
+
+0.4723 0.8815 17
+
+0.4154 0.9096 18
+
+0.3569 0.9341 19
+
+0.2969 0.9549 20
+
+0.7237 0.6901 12
+
+0.6785 0.7346 13
+
+0.8413 0.5406 9
+
+0.8053 0.5929 10
+
+0.766 0.6428 11
+
+0.6306 0.7761 14
+
+0.5801 0.8146 15
+
+0.5272 0.8497 16
+
diff --git a/tests/benchmark_particle_integration_scheme/particles/particles-00001.0000.gnuplot b/tests/benchmark_particle_integration_scheme/particles/particles-00001.0000.gnuplot
new file mode 100644
index 00000000000..caf387eb065
--- /dev/null
+++ b/tests/benchmark_particle_integration_scheme/particles/particles-00001.0000.gnuplot
@@ -0,0 +1,207 @@
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+#
+-0.610517 -0.792008 54
+
+-0.559047 -0.829144 55
+
+-0.505327 -0.862913 56
+
+-0.659479 -0.75169 53
+
+-0.827138 -0.562016 49
+
+-0.789861 -0.6133 50
+
+-0.861178 -0.50838 48
+
+-0.74935 -0.662233 51
+
+-0.705815 -0.708351 52
+
+-0.332867 -0.943026 59
+
+-0.272407 -0.962254 60
+
+-0.449554 -0.893212 57
+
+-0.392009 -0.919997 58
+
+-0.0854358 -0.996339 63
+
+-0.21083 -0.977578 61
+
+-0.0219799 -0.99979 64
+
+-0.148413 -0.988954 62
+
+-0.918531 -0.395287 46
+
+-0.891608 -0.452737 47
+
+-0.961207 -0.275853 44
+
+-0.941769 -0.336272 45
+
+-0.999722 -0.0255667 40
+
+-0.976747 -0.214333 43
+
+-0.988354 -0.151932 42
+
+-0.996071 -0.088929 41
+
+0.229847 -0.973259 68
+
+0.167571 -0.985888 67
+
+0.104715 -0.994487 66
+
+0.041477 -0.99916 65
+
+0.409824 -0.912168 71
+
+0.466848 -0.884319 72
+
+0.291022 -0.956733 69
+
+0.351095 -0.936311 70
+
+0.575032 -0.818172 74
+
+0.625809 -0.780028 75
+
+0.522064 -0.852962 73
+
+0.673933 -0.738744 76
+
+0.719441 -0.69454 77
+
+0.801599 -0.597859 79
+
+0.837904 -0.545757 80
+
+0.762111 -0.647454 78
+
+0.926037 -0.377357 83
+
+0.948191 -0.317862 84
+
+0.966404 -0.257028 85
+
+0.900325 -0.435359 82
+
+0.87082 -0.491544 81
+
+0.980751 -0.195295 86
+
+0.997556 -0.0695921 88
+
+0.99115 -0.132724 87
+
+1.00001 -0.00612132 89
+
+-0.994923 0.101176 38
+
+-0.999351 0.0378753 39
+
+-0.98648 0.164057 37
+
+-0.974067 0.226239 36
+
+-0.937581 0.347827 34
+
+-0.957765 0.287662 35
+
+-0.886039 0.463668 32
+
+-0.913619 0.406535 33
+
+-0.820207 0.572095 30
+
+-0.854807 0.519007 31
+
+-0.782241 0.622932 29
+
+-0.741151 0.671343 28
+
+-0.697062 0.716988 27
+
+-0.548802 0.836001 24
+
+-0.600755 0.799475 25
+
+-0.650195 0.759831 26
+
+-0.380632 0.924732 21
+
+-0.43859 0.898677 22
+
+-0.49468 0.869034 23
+
+-0.321272 0.946985 20
+
+-0.260508 0.965437 19
+
+-0.136234 0.990739 17
+
+-0.198783 0.980014 18
+
+-0.0731287 0.997294 16
+
+-0.00970061 1 15
+
+0.969528 0.245185 93
+
+0.983061 0.183195 92
+
+0.992683 0.120463 91
+
+0.998378 0.0573506 90
+
+0.930573 0.365963 95
+
+0.905562 0.424185 96
+
+0.876867 0.480841 97
+
+0.951981 0.306235 94
+
+0.0537852 0.998553 14
+
+0.179658 0.983726 12
+
+0.116932 0.993146 13
+
+0.241685 0.970335 11
+
+0.362684 0.931938 9
+
+0.302814 0.953077 10
+
+0.477658 0.878529 7
+
+0.420958 0.907042 8
+
+0.727998 0.685656 2
+
+0.682963 0.730436 3
+
+0.808904 0.587952 99
+
+0.77001 0.638061 1
+
+0.844562 0.535492 98
+
+0.808904 0.587952 0
+
+0.585099 0.810991 5
+
+0.532445 0.846524 6
+
+0.635303 0.772195 4
+
diff --git a/tests/benchmark_particle_integration_scheme/screen-output b/tests/benchmark_particle_integration_scheme/screen-output
new file mode 100644
index 00000000000..4219b9fa2db
--- /dev/null
+++ b/tests/benchmark_particle_integration_scheme/screen-output
@@ -0,0 +1,121 @@
+
+Number of active cells: 256 (on 5 levels)
+Number of degrees of freedom: 3,556 (2,178+289+1,089)
+
+*** Timestep 0: t=0 seconds, dt=0 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Writing graphical output: output-benchmark_particle_integration_scheme/solution/solution-00000
+ Writing particle output: output-benchmark_particle_integration_scheme/particles/particles-00000
+
+*** Timestep 1: t=0.00703372 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 2: t=0.0140674 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 3: t=0.0211012 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 4: t=0.0281349 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 5: t=0.0351686 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 6: t=0.0422023 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 7: t=0.049236 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 8: t=0.0562698 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 9: t=0.0633035 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 10: t=0.0703372 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 11: t=0.0773709 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 12: t=0.0844047 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 13: t=0.0914384 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 14: t=0.0984721 seconds, dt=0.00703372 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Number of advected particles: 100
+
+*** Timestep 15: t=0.1 seconds, dt=0.0015279 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+
+ Postprocessing:
+ Writing particle output: output-benchmark_particle_integration_scheme/particles/particles-00001
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/benchmark_particle_integration_scheme/statistics b/tests/benchmark_particle_integration_scheme/statistics
new file mode 100644
index 00000000000..e9258cff717
--- /dev/null
+++ b/tests/benchmark_particle_integration_scheme/statistics
@@ -0,0 +1,26 @@
+# 1: Time step number
+# 2: Time (seconds)
+# 3: Time step size (seconds)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Iterations for temperature solver
+# 8: Visualization file name
+# 9: Number of advected particles
+# 10: Particle file name
+ 0 0.000000000000e+00 0.000000000000e+00 256 2467 1089 0 output-benchmark_particle_integration_scheme/solution/solution-00000 100 output-benchmark_particle_integration_scheme/particles/particles-00000
+ 1 7.033721219977e-03 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+ 2 1.406744243995e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+ 3 2.110116365993e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+ 4 2.813488487991e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+ 5 3.516860609989e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+ 6 4.220232731986e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+ 7 4.923604853984e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+ 8 5.626976975982e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+ 9 6.330349097980e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+10 7.033721219977e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+11 7.737093341975e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+12 8.440465463973e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+13 9.143837585971e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+14 9.847209707968e-02 7.033721219977e-03 256 2467 1089 0 "" 100 ""
+15 1.000000000000e-01 1.527902920317e-03 256 2467 1089 0 "" 100 output-benchmark_particle_integration_scheme/particles/particles-00001