diff --git a/include/aspect/melt.h b/include/aspect/melt.h
index 7f41cd18a8b..3510372816f 100644
--- a/include/aspect/melt.h
+++ b/include/aspect/melt.h
@@ -487,7 +487,7 @@ namespace aspect
        */
       void compute_melt_variables(LinearAlgebra::BlockSparseMatrix &system_matrix,
                                   LinearAlgebra::BlockVector &solution,
-                                  LinearAlgebra::BlockVector &system_rhs);
+                                  LinearAlgebra::BlockVector &system_rhs) const;
 
       /**
        * Return whether this object refers to the porosity field.
diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h
index 1a2dfbcb23c..b8d5fee1ca4 100644
--- a/include/aspect/simulator.h
+++ b/include/aspect/simulator.h
@@ -1444,7 +1444,7 @@ namespace aspect
        * <code>source/simulator/helper_functions.cc</code>.
        */
       void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
-                                            LinearAlgebra::Vector &vec);
+                                            LinearAlgebra::Vector &vec) const;
 
 
       /**
@@ -1476,7 +1476,7 @@ namespace aspect
        * <code>source/simulator/nullspace.cc</code>.
        */
       void remove_nullspace(LinearAlgebra::BlockVector &relevant_dst,
-                            LinearAlgebra::BlockVector &tmp_distributed_stokes);
+                            LinearAlgebra::BlockVector &tmp_distributed_stokes) const;
 
       /**
        * Compute the angular momentum and other rotation properties
@@ -1513,10 +1513,10 @@ namespace aspect
        * This function is implemented in
        * <code>source/simulator/nullspace.cc</code>.
        */
-      void remove_net_angular_momentum( const bool use_constant_density,
-                                        LinearAlgebra::BlockVector &relevant_dst,
-                                        LinearAlgebra::BlockVector &tmp_distributed_stokes,
-                                        const bool limit_to_top_faces = false);
+      void remove_net_angular_momentum(const bool use_constant_density,
+                                       LinearAlgebra::BlockVector &relevant_dst,
+                                       LinearAlgebra::BlockVector &tmp_distributed_stokes,
+                                       const bool limit_to_top_faces = false) const;
 
       /**
        * Offset the boundary id of all faces located on an outflow boundary
@@ -1549,9 +1549,9 @@ namespace aspect
        * This function is implemented in
        * <code>source/simulator/nullspace.cc</code>.
        */
-      void remove_net_linear_momentum( const bool use_constant_density,
-                                       LinearAlgebra::BlockVector &relevant_dst,
-                                       LinearAlgebra::BlockVector &tmp_distributed_stokes);
+      void remove_net_linear_momentum(const bool use_constant_density,
+                                      LinearAlgebra::BlockVector &relevant_dst,
+                                      LinearAlgebra::BlockVector &tmp_distributed_stokes) const;
 
       /**
        * Compute the maximal velocity throughout the domain. This is needed to
diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index bbb700565e1..6bd4a311d4d 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -823,7 +823,7 @@ namespace aspect
 
   template <int dim>
   void Simulator<dim>::interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
-                                                        LinearAlgebra::Vector &vec)
+                                                        LinearAlgebra::Vector &vec) const
   {
     Assert(introspection.block_indices.velocities == 0, ExcNotImplemented());
 
@@ -2771,7 +2771,7 @@ namespace aspect
   template double Simulator<dim>::compute_initial_stokes_residual(); \
   template bool Simulator<dim>::stokes_matrix_depends_on_solution() const; \
   template bool Simulator<dim>::stokes_A_block_is_symmetric() const; \
-  template void Simulator<dim>::interpolate_onto_velocity_system(const TensorFunction<1,dim> &func, LinearAlgebra::Vector &vec);\
+  template void Simulator<dim>::interpolate_onto_velocity_system(const TensorFunction<1,dim> &func, LinearAlgebra::Vector &vec) const;\
   template void Simulator<dim>::apply_limiter_to_dg_solutions(const AdvectionField &advection_field); \
   template void Simulator<dim>::compute_reactions(); \
   template void Simulator<dim>::initialize_current_linearization_point (); \
diff --git a/source/simulator/melt.cc b/source/simulator/melt.cc
index c5110d07924..78e9224a76b 100644
--- a/source/simulator/melt.cc
+++ b/source/simulator/melt.cc
@@ -1074,7 +1074,7 @@ namespace aspect
   MeltHandler<dim>::
   compute_melt_variables(LinearAlgebra::BlockSparseMatrix &system_matrix,
                          LinearAlgebra::BlockVector &solution,
-                         LinearAlgebra::BlockVector &system_rhs)
+                         LinearAlgebra::BlockVector &system_rhs) const
   {
     if (!this->include_melt_transport())
       return;
diff --git a/source/simulator/nullspace.cc b/source/simulator/nullspace.cc
index e9644fdab2b..92e696eb1f2 100644
--- a/source/simulator/nullspace.cc
+++ b/source/simulator/nullspace.cc
@@ -207,7 +207,7 @@ namespace aspect
 
   template <int dim>
   void Simulator<dim>::remove_nullspace(LinearAlgebra::BlockVector &relevant_dst,
-                                        LinearAlgebra::BlockVector &tmp_distributed_stokes)
+                                        LinearAlgebra::BlockVector &tmp_distributed_stokes) const
   {
     if (parameters.nullspace_removal & NullspaceRemoval::angular_momentum)
       {
@@ -244,9 +244,9 @@ namespace aspect
   }
 
   template <int dim>
-  void Simulator<dim>::remove_net_linear_momentum( const bool use_constant_density,
-                                                   LinearAlgebra::BlockVector &relevant_dst,
-                                                   LinearAlgebra::BlockVector &tmp_distributed_stokes )
+  void Simulator<dim>::remove_net_linear_momentum(const bool use_constant_density,
+                                                  LinearAlgebra::BlockVector &relevant_dst,
+                                                  LinearAlgebra::BlockVector &tmp_distributed_stokes) const
   {
     Assert(introspection.block_indices.velocities != introspection.block_indices.pressure,
            ExcNotImplemented());
@@ -467,10 +467,10 @@ namespace aspect
 
 
   template <int dim>
-  void Simulator<dim>::remove_net_angular_momentum( const bool use_constant_density,
-                                                    LinearAlgebra::BlockVector &relevant_dst,
-                                                    LinearAlgebra::BlockVector &tmp_distributed_stokes,
-                                                    const bool limit_to_top_faces)
+  void Simulator<dim>::remove_net_angular_momentum(const bool use_constant_density,
+                                                   LinearAlgebra::BlockVector &relevant_dst,
+                                                   LinearAlgebra::BlockVector &tmp_distributed_stokes,
+                                                   const bool limit_to_top_faces) const
   {
     Assert(introspection.block_indices.velocities != introspection.block_indices.pressure,
            ExcNotImplemented());
@@ -513,7 +513,7 @@ namespace aspect
 {
 #define INSTANTIATE(dim) \
   template struct RotationProperties<dim>; \
-  template void Simulator<dim>::remove_nullspace (LinearAlgebra::BlockVector &,LinearAlgebra::BlockVector &vector); \
+  template void Simulator<dim>::remove_nullspace (LinearAlgebra::BlockVector &,LinearAlgebra::BlockVector &vector) const; \
   template void Simulator<dim>::setup_nullspace_constraints (AffineConstraints<double> &);
 
   ASPECT_INSTANTIATE(INSTANTIATE)