diff --git a/include/io/vtk_writer.h b/include/io/vtk_writer.h
index eb0591d7e..549114d78 100644
--- a/include/io/vtk_writer.h
+++ b/include/io/vtk_writer.h
@@ -74,10 +74,11 @@ class VtkWriter {
   //! \param[in] mpi_size Number of MPI tasks
   //! \param[in] step Current time step
   //! \param[in] max_steps Maximum number of steps in the simulation
+  //! \param[in] step_size Physical time step size
   //! \param[in] ncomponents Number of components to write
   void write_parallel_vtk(const std::string& filename,
                           const std::string& attribute, int mpi_size,
-                          unsigned step, unsigned max_steps,
+                          unsigned step, unsigned max_steps, double step_size,
                           unsigned ncomponents = 3);
 
  private:
diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc
index f829b5f95..9df2b65cb 100644
--- a/include/solvers/mpm_base.tcc
+++ b/include/solvers/mpm_base.tcc
@@ -541,7 +541,7 @@ void mpm::MPMBase<Tdim>::write_vtk(mpm::Index step, mpm::Index max_steps) {
                                .string();
 
       vtk_writer->write_parallel_vtk(parallel_file, attribute, mpi_size, step,
-                                     max_steps);
+                                     max_steps, dt_);
     }
 #endif
   }
@@ -562,7 +562,7 @@ void mpm::MPMBase<Tdim>::write_vtk(mpm::Index step, mpm::Index max_steps) {
                                .string();
 
       vtk_writer->write_parallel_vtk(parallel_file, attribute, mpi_size, step,
-                                     max_steps);
+                                     max_steps, dt_);
     }
 #endif
   }
@@ -583,7 +583,7 @@ void mpm::MPMBase<Tdim>::write_vtk(mpm::Index step, mpm::Index max_steps) {
                                .string();
 
       vtk_writer->write_parallel_vtk(parallel_file, attribute, mpi_size, step,
-                                     max_steps, 9);
+                                     max_steps, dt_, 9);
     }
 #endif
   }
@@ -609,7 +609,7 @@ void mpm::MPMBase<Tdim>::write_vtk(mpm::Index step, mpm::Index max_steps) {
                                  .string();
         unsigned ncomponents = 1;
         vtk_writer->write_parallel_vtk(parallel_file, phase_attribute, mpi_size,
-                                       step, max_steps, ncomponents);
+                                       step, max_steps, dt_, ncomponents);
       }
 #endif
     }
diff --git a/src/io/vtk_writer.cc b/src/io/vtk_writer.cc
index 6a6a51e9b..1239ddc1b 100644
--- a/src/io/vtk_writer.cc
+++ b/src/io/vtk_writer.cc
@@ -1,4 +1,5 @@
 #include "vtk_writer.h"
+#include <boost/filesystem.hpp>
 
 #ifdef USE_VTK
 
@@ -223,7 +224,7 @@ void VtkWriter::write_mesh(
 void VtkWriter::write_parallel_vtk(const std::string& filename,
                                    const std::string& attribute, int mpi_size,
                                    unsigned step, unsigned max_steps,
-                                   unsigned ncomponents) {
+                                   double step_size, unsigned ncomponents) {
 
   // If the number of components is 1, set as scalar or vector / tensor
   std::string data_type;
@@ -278,6 +279,46 @@ void VtkWriter::write_parallel_vtk(const std::string& filename,
   pvtk.open(filename);
   pvtk << ppolydata;
   pvtk.close();
+
+  // Write parallel grouping VTK file
+
+  std::string output_path = filename.substr(0, filename.find_last_of("\\/"));
+  std::ofstream group_vtk;
+  std::string group_filename = output_path + "/" + attribute + ".pvd";
+  std::string group_data;
+
+  std::stringstream group_parts_file;
+  group_parts_file.str(std::string());
+  group_parts_file << attribute;
+  group_parts_file.fill('0');
+  int digits = log10(max_steps) + 1;
+  group_parts_file.width(digits);
+  group_parts_file << step;
+  group_parts_file << ".pvtp";
+
+  boost::filesystem::path file_check(group_filename);
+
+  if (boost::filesystem::exists(file_check)) {
+    group_vtk.open(group_filename, std::fstream::app);
+    group_data = "\t<DataSet timestep=\"" + std::to_string(step * step_size) +
+                 "\" file=\"./" + group_parts_file.str() + "\"/>\n";
+    group_vtk << group_data;
+  } else {
+    group_vtk.open(group_filename);
+    group_data =
+        "<?xml version=\"1.0\"?>\n<VTKFile type=\"Collection\" version=\"0.1\" "
+        "byte_order=\"LittleEndian\">\n<Collection>\n\t<DataSet timestep=\"" +
+        std::to_string(step * step_size) + "\" file=\"./" +
+        group_parts_file.str() + "\"/>\n";
+    group_vtk << group_data;
+  }
+
+  if (step == max_steps - 1) {
+    std::string closing = "</Collection>\n</VTKFile>";
+    group_vtk << closing;
+  }
+
+  group_vtk.close();
 }
 
 #endif
diff --git a/tests/io/vtk_writer_test.cc b/tests/io/vtk_writer_test.cc
index da07b83cb..a0f6c1ccd 100644
--- a/tests/io/vtk_writer_test.cc
+++ b/tests/io/vtk_writer_test.cc
@@ -51,8 +51,9 @@ TEST_CASE("VTK Writer is checked", "[vtk][writer]") {
     int mpi_size = 2;
     unsigned step = 1000;
     unsigned max_steps = 10000;
+    double step_size = 0.5;
     vtk_writer->write_parallel_vtk(parallel_vtk_file, attribute, mpi_size, step,
-                                   max_steps);
+                                   max_steps, step_size);
 
     // Check file data
     std::string ppolydata =
@@ -98,8 +99,9 @@ TEST_CASE("VTK Writer is checked", "[vtk][writer]") {
     int mpi_size = 2;
     unsigned step = 1000;
     unsigned max_steps = 10000;
+    double step_size = 0.5;
     vtk_writer->write_parallel_vtk(parallel_vtk_file, attribute, mpi_size, step,
-                                   max_steps, 1);
+                                   max_steps, step_size, 1);
 
     // Check file data
     std::string ppolydata =