diff --git a/src/state/ageinglawstateupdater.cc b/src/state/ageinglawstateupdater.cc
index ec69fac30fa339520e8ec6b6c719860d43c22399..51a916267a30859d64cd477c8e1744a9dc2f5b79 100644
--- a/src/state/ageinglawstateupdater.cc
+++ b/src/state/ageinglawstateupdater.cc
@@ -46,3 +46,9 @@ void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
     ScalarVector &_alpha) {
   _alpha = alpha;
 }
+
+template <class ScalarVector, class Vector>
+std::shared_ptr<StateUpdater<ScalarVector, Vector>>
+AgeingLawStateUpdater<ScalarVector, Vector>::clone() const {
+  return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(*this);
+}
diff --git a/src/state/ageinglawstateupdater.hh b/src/state/ageinglawstateupdater.hh
index c8ae85c7e348a4598d340f93f731289c386be328..67ebab25ee7548916ecaddb6c1de3799f3f1604f 100644
--- a/src/state/ageinglawstateupdater.hh
+++ b/src/state/ageinglawstateupdater.hh
@@ -16,6 +16,8 @@ class AgeingLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
   void solve(Vector const &velocity_field) override;
   void extractAlpha(ScalarVector &) override;
 
+  std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const;
+
 private:
   ScalarVector alpha_o;
   ScalarVector alpha;
diff --git a/src/state/sliplawstateupdater.cc b/src/state/sliplawstateupdater.cc
index 2b3e643c39051425959e014e598666a07118700b..acffc6f217fbf2cb8bbd1cea7e16b41ff58acb86 100644
--- a/src/state/sliplawstateupdater.cc
+++ b/src/state/sliplawstateupdater.cc
@@ -35,3 +35,9 @@ void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha(
     ScalarVector &_alpha) {
   _alpha = alpha;
 }
+
+template <class ScalarVector, class Vector>
+std::shared_ptr<StateUpdater<ScalarVector, Vector>>
+SlipLawStateUpdater<ScalarVector, Vector>::clone() const {
+  return std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(*this);
+}
diff --git a/src/state/sliplawstateupdater.hh b/src/state/sliplawstateupdater.hh
index e12e004ca266f9cb386f34033f46b8930abe41d6..28838aec28f128f97770734b1b2fceff9464b927 100644
--- a/src/state/sliplawstateupdater.hh
+++ b/src/state/sliplawstateupdater.hh
@@ -16,6 +16,8 @@ class SlipLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
   void solve(Vector const &velocity_field) override;
   void extractAlpha(ScalarVector &) override;
 
+  std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const;
+
 private:
   ScalarVector alpha_o;
   ScalarVector alpha;
diff --git a/src/state/stateupdater.hh b/src/state/stateupdater.hh
index 697662d32e8bf877441fed9085c7ccaa0bba087c..9e4495af34e2444345f1f718dc4e077bba2e0bd2 100644
--- a/src/state/stateupdater.hh
+++ b/src/state/stateupdater.hh
@@ -9,6 +9,8 @@ template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
   void virtual setup(double _tau) = 0;
   void virtual solve(Vector const &velocity_field) = 0;
   void virtual extractAlpha(ScalarVector &alpha) = 0;
+
+  std::shared_ptr<StateUpdater<ScalarVector, Vector>> virtual clone() const = 0;
 };
 
 #endif
diff --git a/src/timestepping.hh b/src/timestepping.hh
index 71d99b69de6ee66fa8c2aae9f8299737941af5ba..1221c585f6c5a4b728d3ade1f3725642f9536ca3 100644
--- a/src/timestepping.hh
+++ b/src/timestepping.hh
@@ -20,6 +20,10 @@ class TimeSteppingScheme {
   void virtual extractVelocity(Vector &velocity) const = 0;
   void virtual extractRelativeVelocity(Vector &velocity) const = 0;
   void virtual extractOldVelocity(Vector &velocity) const = 0;
+
+  std::shared_ptr<
+      TimeSteppingScheme<Vector, Matrix, Function, dim>> virtual clone()
+      const = 0;
 };
 
 #include "timestepping/newmark.hh"
diff --git a/src/timestepping/backward_euler.cc b/src/timestepping/backward_euler.cc
index 57f24bd3998397eea314d48d775748d094d96a79..a39f0e5415a0ec0faf028badcd489a1444e51477 100644
--- a/src/timestepping/backward_euler.cc
+++ b/src/timestepping/backward_euler.cc
@@ -147,3 +147,9 @@ void BackwardEuler<Vector, Matrix, Function, dim>::extractOldVelocity(
     Vector &velocity) const {
   velocity = v_o;
 }
+
+template <class Vector, class Matrix, class Function, size_t dim>
+std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dim>>
+BackwardEuler<Vector, Matrix, Function, dim>::clone() const {
+  return std::make_shared<BackwardEuler<Vector, Matrix, Function, dim>>(*this);
+}
diff --git a/src/timestepping/backward_euler.hh b/src/timestepping/backward_euler.hh
index e6291fd69e9d08a50c2d518cc0d517faea87010e..510a149f8a84b9b52dba82a5fee28c7a0a6f17b9 100644
--- a/src/timestepping/backward_euler.hh
+++ b/src/timestepping/backward_euler.hh
@@ -21,6 +21,9 @@ class BackwardEuler : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
   void extractRelativeVelocity(Vector &) const override;
   void extractOldVelocity(Vector &) const override;
 
+  std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dim>> clone()
+      const;
+
 private:
   Matrix const &A;
   Matrix const &M;
diff --git a/src/timestepping/newmark.cc b/src/timestepping/newmark.cc
index 518a69ce1ab53c090dffd78537ca0c2df223c11a..51eeb7a36f79283e4fd9ed928203b976499c7bf4 100644
--- a/src/timestepping/newmark.cc
+++ b/src/timestepping/newmark.cc
@@ -166,3 +166,9 @@ void Newmark<Vector, Matrix, Function, dim>::extractOldVelocity(
     Vector &velocity) const {
   velocity = v_o;
 }
+
+template <class Vector, class Matrix, class Function, size_t dim>
+std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dim>>
+Newmark<Vector, Matrix, Function, dim>::clone() const {
+  return std::make_shared<Newmark<Vector, Matrix, Function, dim>>(*this);
+}
diff --git a/src/timestepping/newmark.hh b/src/timestepping/newmark.hh
index bfae966491ab556ea6df287f846ea7cb2c0612ad..0d711447b5b94c33f5d0299ef7c7529a8564a1e8 100644
--- a/src/timestepping/newmark.hh
+++ b/src/timestepping/newmark.hh
@@ -22,6 +22,9 @@ class Newmark : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
   void extractRelativeVelocity(Vector &) const override;
   void extractOldVelocity(Vector &) const override;
 
+  std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dim>> clone()
+      const;
+
 private:
   Matrix const &A;
   Matrix const &M;