diff --git a/include/aspect/compat.h b/include/aspect/compat.h
index ead5c2e4cb9..40abd6e169f 100644
--- a/include/aspect/compat.h
+++ b/include/aspect/compat.h
@@ -444,6 +444,15 @@ namespace aspect
                                const double      outer_radius);
 #endif
 
+  // deal.II 9.8 made ReferenceCell a template class, whereas older versions
+  // had it as a non-template class. This is a problem.
+  // Rather than litter our own code base with #ifdefs, we can just define the
+  // templated class variant here for older deal.II versions, and then we can
+  // use the same code in all versions.
+#if !DEAL_II_VERSION_GTE(9,8,0)
+  template <int dim> using ReferenceCell = dealii::ReferenceCell;
+#endif
+
 }
 
 #endif
diff --git a/include/aspect/simulator/assemblers/interface.h b/include/aspect/simulator/assemblers/interface.h
index 5ca60b0152f..a6aa3a83f82 100644
--- a/include/aspect/simulator/assemblers/interface.h
+++ b/include/aspect/simulator/assemblers/interface.h
@@ -487,7 +487,6 @@ namespace aspect
    */
   namespace Assemblers
   {
-
     /**
      * For a reference cell (which is typically obtained by asking the finite
      * element to be used), determine how many interface matrices are needed.
@@ -497,16 +496,18 @@ namespace aspect
      * accommodates the fact that the neighbors of a cell can all be refined,
      * though they can only be refined once.
      */
+    template <int dim>
     unsigned int
-    n_interface_matrices (const ReferenceCell &reference_cell);
+    n_interface_matrices (const ReferenceCell<dim> &reference_cell);
 
     /**
      * For a given reference cell, and a given face we are currently
      * assembling on, return which element of an array of size
      * `n_interface_matrices(reference_cell)` to use.
      */
+    template <int dim>
     unsigned int
-    nth_interface_matrix (const ReferenceCell &reference_cell,
+    nth_interface_matrix (const ReferenceCell<dim> &reference_cell,
                           const unsigned int face);
 
     /**
@@ -514,8 +515,9 @@ namespace aspect
      * currently assembling on, return which element of an array of size
      * `n_interface_matrices(reference_cell)` to use.
      */
+    template <int dim>
     unsigned int
-    nth_interface_matrix (const ReferenceCell &reference_cell,
+    nth_interface_matrix (const ReferenceCell<dim> &reference_cell,
                           const unsigned int face,
                           const unsigned int sub_face);
 
diff --git a/source/simulator/assemblers/interface.cc b/source/simulator/assemblers/interface.cc
index c5581d40a7a..ff29b582269 100644
--- a/source/simulator/assemblers/interface.cc
+++ b/source/simulator/assemblers/interface.cc
@@ -501,8 +501,9 @@ namespace aspect
 
   namespace Assemblers
   {
+    template <int dim>
     unsigned int
-    n_interface_matrices (const ReferenceCell &reference_cell)
+    n_interface_matrices (const ReferenceCell<dim> &reference_cell)
     {
       // The current implementation assumes that all faces are
       // the same; so no wedges or pyramids please.
@@ -520,8 +521,9 @@ namespace aspect
 
 
 
+    template <int dim>
     unsigned int
-    nth_interface_matrix (const ReferenceCell &reference_cell,
+    nth_interface_matrix (const ReferenceCell<dim> &reference_cell,
                           const unsigned int face)
     {
       AssertIndexRange (face, reference_cell.n_faces());
@@ -531,8 +533,9 @@ namespace aspect
 
 
 
+    template <int dim>
     unsigned int
-    nth_interface_matrix (const ReferenceCell &reference_cell,
+    nth_interface_matrix (const ReferenceCell<dim> &reference_cell,
                           const unsigned int face,
                           const unsigned int sub_face)
     {
@@ -649,7 +652,23 @@ namespace aspect
     template class Interface<dim>; \
     template class AdvectionStabilizationInterface<dim>; \
     template class Manager<dim>; \
+    \
+    template \
+    unsigned int \
+    n_interface_matrices (const ReferenceCell<dim> &reference_cell); \
+    \
+    template \
+    unsigned int \
+    nth_interface_matrix (const ReferenceCell<dim> &reference_cell, \
+                          const unsigned int face); \
+    \
+    template \
+    unsigned int \
+    nth_interface_matrix (const ReferenceCell<dim> &reference_cell, \
+                          const unsigned int face, \
+                          const unsigned int sub_face); \
   }
+
   ASPECT_INSTANTIATE(INSTANTIATE)
 
 #undef INSTANTIATE
diff --git a/source/simulator/introspection.cc b/source/simulator/introspection.cc
index 495539139bc..4824ed8e0d8 100644
--- a/source/simulator/introspection.cc
+++ b/source/simulator/introspection.cc
@@ -140,24 +140,33 @@ namespace aspect
     template <int dim>
     typename Introspection<dim>::Quadratures
     setup_quadratures (const Parameters<dim> &parameters,
-                       const ReferenceCell reference_cell)
+                       const ReferenceCell<dim> reference_cell)
     {
+      const auto get_gauss_type_quadrature = [&reference_cell](const unsigned int degree)
+      {
+#if DEAL_II_VERSION_GTE(9,8,0)
+        return reference_cell.get_gauss_type_quadrature(degree);
+#else
+        return reference_cell.get_gauss_type_quadrature<dim>(degree);
+#endif
+      };
+
       typename Introspection<dim>::Quadratures quadratures;
 
-      quadratures.velocities = reference_cell.get_gauss_type_quadrature<dim>(parameters.stokes_velocity_degree+1);
-      quadratures.pressure = reference_cell.get_gauss_type_quadrature<dim>(parameters.use_equal_order_interpolation_for_stokes
-                                                                           ?
-                                                                           parameters.stokes_velocity_degree+1
-                                                                           :
-                                                                           parameters.stokes_velocity_degree);
-      quadratures.temperature = reference_cell.get_gauss_type_quadrature<dim>(parameters.temperature_degree+1);
-      quadratures.compositional_field_max = reference_cell.get_gauss_type_quadrature<dim>(parameters.max_composition_degree+1);
+      quadratures.velocities = get_gauss_type_quadrature(parameters.stokes_velocity_degree+1);
+      quadratures.pressure = get_gauss_type_quadrature(parameters.use_equal_order_interpolation_for_stokes
+                                                       ?
+                                                       parameters.stokes_velocity_degree+1
+                                                       :
+                                                       parameters.stokes_velocity_degree);
+      quadratures.temperature = get_gauss_type_quadrature(parameters.temperature_degree+1);
+      quadratures.compositional_field_max = get_gauss_type_quadrature(parameters.max_composition_degree+1);
       for (const auto &degree: parameters.composition_degrees)
-        quadratures.compositional_fields.emplace_back(reference_cell.get_gauss_type_quadrature<dim>(degree+1));
-      quadratures.system = reference_cell.get_gauss_type_quadrature<dim>(std::max({parameters.stokes_velocity_degree,
-                                                                                   parameters.temperature_degree,
-                                                                                   parameters.max_composition_degree
-                                                                                  }) + 1);
+        quadratures.compositional_fields.emplace_back(get_gauss_type_quadrature(degree+1));
+      quadratures.system = get_gauss_type_quadrature(std::max({parameters.stokes_velocity_degree,
+                                                               parameters.temperature_degree,
+                                                               parameters.max_composition_degree
+                                                              }) + 1);
 
       return quadratures;
     }
@@ -167,19 +176,29 @@ namespace aspect
     template <int dim>
     typename Introspection<dim>::FaceQuadratures
     setup_face_quadratures (const Parameters<dim> &parameters,
-                            const ReferenceCell reference_cell)
+                            const ReferenceCell<dim> reference_cell)
     {
+      const auto get_gauss_type_face_quadrature = [&reference_cell](const unsigned int degree)
+      {
+#if DEAL_II_VERSION_GTE(9,8,0)
+        return reference_cell.face_reference_cell(0).get_gauss_type_quadrature(degree);
+#else
+        return reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(degree);
+#endif
+      };
+
+
       typename Introspection<dim>::FaceQuadratures quadratures;
 
-      quadratures.velocities = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(parameters.stokes_velocity_degree+1);
-      quadratures.pressure = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(parameters.use_equal_order_interpolation_for_stokes
-                             ?
-                             parameters.stokes_velocity_degree+1
-                             :
-                             parameters.stokes_velocity_degree);
-      quadratures.temperature = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(parameters.temperature_degree+1);
-      quadratures.compositional_fields = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(parameters.max_composition_degree+1);
-      quadratures.system = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(
+      quadratures.velocities = get_gauss_type_face_quadrature(parameters.stokes_velocity_degree+1);
+      quadratures.pressure = get_gauss_type_face_quadrature(parameters.use_equal_order_interpolation_for_stokes
+                                                            ?
+                                                            parameters.stokes_velocity_degree+1
+                                                            :
+                                                            parameters.stokes_velocity_degree);
+      quadratures.temperature = get_gauss_type_face_quadrature(parameters.temperature_degree+1);
+      quadratures.compositional_fields = get_gauss_type_face_quadrature(parameters.max_composition_degree+1);
+      quadratures.system = get_gauss_type_face_quadrature(
                              std::max({parameters.stokes_velocity_degree, parameters.temperature_degree, parameters.max_composition_degree}) + 1);
 
       return quadratures;