diff --git a/.gitignore b/.gitignore
index 935ed9c3483a235a260a05822ea82f628e16f958..427847af023de7e0b4e25b93e408b39744916e9f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -24,3 +24,5 @@ Tests/Data/Parabolic/HT/InvalidProjectFiles/*.prj
 web/.netlify
 
 CMakeUserPresets.json
+
+.snakemake
diff --git a/Applications/Utils/Tests.cmake b/Applications/Utils/Tests.cmake
index 6c5ab390489f993fa2cdc5a40a3f159f0743e0b3..a1ac37c19da430f3a2922a6909c03def3687d3b3 100644
--- a/Applications/Utils/Tests.cmake
+++ b/Applications/Utils/Tests.cmake
@@ -301,14 +301,17 @@ AddTest(
 
 if(SNAKEMAKE AND NOT OGS_USE_MPI)
     add_test(NAME snakemake_ExtractBoundary
-        COMMAND ${SNAKEMAKE}
-                 -j 1
-                 --configfile ${PROJECT_BINARY_DIR}/buildinfo.yaml
-                 -s ${CMAKE_CURRENT_SOURCE_DIR}/ExtractBoundary.smk
+        COMMAND ${SNAKEMAKE} -j 1
+            --configfile ${PROJECT_BINARY_DIR}/buildinfo.yaml
+            -s ${CMAKE_CURRENT_SOURCE_DIR}/ExtractBoundary.smk
     )
-endif()
-if(SNAKEMAKE)
-    add_dependencies(ctest ExtractBoundary)
+
+    add_test(NAME snakemake_VoxelGridFromLayers
+        COMMAND ${SNAKEMAKE} -j 1
+            --configfile ${PROJECT_BINARY_DIR}/buildinfo.yaml
+            -s ${CMAKE_CURRENT_SOURCE_DIR}/VoxelGridFromLayers.smk
+    )
+    add_dependencies(ctest ExtractBoundary Layers2Grid AddFaultToVoxelGrid)
 endif()
 
 AddTest(
@@ -545,33 +548,3 @@ AddTest(
     DIFF_DATA
     PrismBHE_elev.vtu
 )
-
-MeshTest(
-    NAME Layers2Grid_Iso_Test
-    PATH MeshLib
-    WORKING_DIRECTORY ${Data_SOURCE_DIR}/MeshLib
-    EXECUTABLE Layers2Grid
-    EXECUTABLE_ARGS -i AREHS_test_layers.txt -o ${Data_BINARY_DIR}/MeshLib/AREHS_test_iso.vtu -x 500
-    REQUIREMENTS NOT OGS_USE_MPI
-    DIFF_DATA AREHS_test_iso.vtu AREHS_test_iso.vtu 1e-16
-)
-
-MeshTest(
-    NAME Layers2Grid_Test
-    PATH MeshLib
-    WORKING_DIRECTORY ${Data_SOURCE_DIR}/MeshLib
-    EXECUTABLE Layers2Grid
-    EXECUTABLE_ARGS -i AREHS_test_layers.txt -o ${Data_BINARY_DIR}/MeshLib/AREHS_test.vtu -x 500 -y 300 -z 100
-    REQUIREMENTS NOT OGS_USE_MPI
-    DIFF_DATA AREHS_test.vtu AREHS_test.vtu 1e-16
-)
-
-MeshTest(
-    NAME AddFaultToVoxelGrid_Test
-    PATH MeshLib
-    WORKING_DIRECTORY ${Data_SOURCE_DIR}/MeshLib
-    EXECUTABLE AddFaultToVoxelGrid
-    EXECUTABLE_ARGS -i AREHS_test.vtu -f AREHS_fault.vtu -o ${Data_BINARY_DIR}/MeshLib/AREHS_test_fault.vtu
-    REQUIREMENTS NOT OGS_USE_MPI
-    DIFF_DATA AREHS_test_fault.vtu AREHS_test_fault.vtu 1e-16
-)
diff --git a/Applications/Utils/VoxelGridFromLayers.smk b/Applications/Utils/VoxelGridFromLayers.smk
new file mode 100644
index 0000000000000000000000000000000000000000..921c9e153094dd9de071961193e4c80525ed2c76
--- /dev/null
+++ b/Applications/Utils/VoxelGridFromLayers.smk
@@ -0,0 +1,77 @@
+# Usage, e.g.:
+#   snakemake -s VoxelGridFromLayers.smk -j 1 --configfile $HOME/code/ogs6/build/buildinfo.yaml
+#
+# buildinfo.yaml contains variables such as Data_BINARY_DIR
+
+import os
+os.environ["PATH"] += os.pathsep + os.pathsep.join([config['BIN_DIR']])
+workdir: f"{config['Data_SOURCE_DIR']}/MeshLib"
+out_dir = f"{config['Data_BINARY_DIR']}/MeshLib"
+
+rule all:
+    input:
+        f"{out_dir}/AREHS_test_diff_geometry.out",
+        f"{out_dir}/AREHS_test_fault_diff_material_ids.out",
+        f"{out_dir}/AREHS_test_fault_diff_geometry.out",
+        f"{out_dir}/AREHS_test_iso_diff_geometry.out"
+
+rule layers_to_grid:
+    input:
+        "AREHS_test_layers.txt"
+    output:
+        f"{out_dir}/AREHS_test.vtu"
+    shell:
+        "Layers2Grid -i {input} -o {output} -x 500 -y 300 -z 100"
+
+rule vtkdiff_grid_geometry:
+    input:
+        out = rules.layers_to_grid.output,
+        ref = "AREHS_test.vtu"
+    output:
+        f"{out_dir}/AREHS_test_diff_geometry.out"
+    shell:
+        "vtkdiff -m {input.out} {input.ref} > {output}"
+
+rule add_fault_to_grid:
+    input:
+        grid = rules.layers_to_grid.output,
+        fault = "AREHS_fault.vtu"
+    output:
+        f"{out_dir}/AREHS_test_fault.vtu"
+    shell:
+        "AddFaultToVoxelGrid -i {input.grid} -f {input.fault} -o {output}"
+
+rule vtkdiff_material_ids:
+    input:
+        out = rules.add_fault_to_grid.output,
+        ref = "AREHS_test_fault.vtu"
+    output:
+        f"{out_dir}/AREHS_test_fault_diff_material_ids.out"
+    shell:
+        "vtkdiff -a MaterialIDs -b MaterialIDs {input.out} {input.ref} > {output}"
+
+rule vtkdiff_fault_geometry:
+    input:
+        out = rules.add_fault_to_grid.output,
+        ref = "AREHS_test_fault.vtu"
+    output:
+        f"{out_dir}/AREHS_test_fault_diff_geometry.out"
+    shell:
+        "vtkdiff -m {input.out} {input.ref} > {output}"
+
+rule layers_to_grid_iso:
+    input:
+        "AREHS_test_layers.txt"
+    output:
+        f"{out_dir}/AREHS_test_iso.vtu"
+    shell:
+        "Layers2Grid -i {input} -o {output} -x 500"
+
+rule vtkdiff_grid_iso_geometry:
+    input:
+        out = rules.layers_to_grid_iso.output,
+        ref = "AREHS_test_iso.vtu"
+    output:
+        f"{out_dir}/AREHS_test_iso_diff_geometry.out"
+    shell:
+        "vtkdiff -m {input.out} {input.ref} > {output}"