diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h
index 4d5c99386d4..ece36e11799 100644
--- a/include/aspect/simulator.h
+++ b/include/aspect/simulator.h
@@ -1545,7 +1545,9 @@ namespace aspect
        * This function is implemented in
        * <code>source/simulator/helper_functions.cc</code>.
        */
-      void replace_outflow_boundary_ids(const unsigned int boundary_id_offset);
+      void replace_outflow_boundary_ids(const unsigned int boundary_id_offset,
+                                        const bool is_composition,
+                                        const unsigned int composition_index);
 
       /**
        * Undo the offset of the boundary ids done in replace_outflow_boundary_ids
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index e8173dae562..81d4bd85fc5 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -695,7 +695,7 @@ namespace aspect
     // so we want to offset them by 128 and not allow more than 128 boundary ids.
     const unsigned int boundary_id_offset = 128;
     if (!boundary_temperature_manager.allows_fixed_temperature_on_outflow_boundaries())
-      replace_outflow_boundary_ids(boundary_id_offset);
+      replace_outflow_boundary_ids(boundary_id_offset, false, numbers::invalid_unsigned_int);
 
     // if using continuous temperature FE, do the same for the temperature variable:
     // evaluate the current boundary temperature and add these constraints as well
@@ -730,41 +730,43 @@ namespace aspect
     // update the composition boundary condition.
     boundary_composition_manager.update();
 
-    // If we do not want to prescribe Dirichlet boundary conditions on outflow boundaries,
-    // use the same trick for marking up outflow boundary conditions for compositional fields
-    // as we did above already for the temperature.
-    if (!boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries())
-      replace_outflow_boundary_ids(boundary_id_offset);
-
     // now do the same for the composition variables:
     {
       // obtain the boundary indicators that belong to Dirichlet-type
       // composition boundary conditions and interpolate the composition
       // there
       for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
-        if (parameters.use_discontinuous_composition_discretization[c] == false)
-          for (const auto p : boundary_composition_manager.get_fixed_composition_boundary_indicators())
-            {
-              VectorFunctionFromScalarFunctionObject<dim> vector_function_object(
-                [&] (const Point<dim> &x) -> double
+        {
+          // If we do not want to prescribe Dirichlet boundary conditions on outflow boundaries,
+          // use the same trick for marking up outflow boundary conditions for compositional fields
+          // as we did above already for the temperature.
+          if (!boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries())
+            replace_outflow_boundary_ids(boundary_id_offset, true, c);
+
+          if (parameters.use_discontinuous_composition_discretization[c] == false)
+            for (const auto p : boundary_composition_manager.get_fixed_composition_boundary_indicators())
               {
-                return boundary_composition_manager.boundary_composition(p, x, c);
-              },
-              introspection.component_masks.compositional_fields[c].first_selected_component(),
-              introspection.n_components);
-
-              VectorTools::interpolate_boundary_values (*mapping,
-                                                        dof_handler,
-                                                        p,
-                                                        vector_function_object,
-                                                        new_current_constraints,
-                                                        introspection.component_masks.compositional_fields[c]);
-            }
+                VectorFunctionFromScalarFunctionObject<dim> vector_function_object(
+                  [&] (const Point<dim> &x) -> double
+                {
+                  return boundary_composition_manager.boundary_composition(p, x, c);
+                },
+                introspection.component_masks.compositional_fields[c].first_selected_component(),
+                introspection.n_components);
+
+                VectorTools::interpolate_boundary_values (*mapping,
+                                                          dof_handler,
+                                                          p,
+                                                          vector_function_object,
+                                                          new_current_constraints,
+                                                          introspection.component_masks.compositional_fields[c]);
+
+                if (!boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries())
+                  restore_outflow_boundary_ids(boundary_id_offset);
+              }
+        }
     }
 
-    if (!boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries())
-      restore_outflow_boundary_ids(boundary_id_offset);
-
     if (parameters.include_melt_transport)
       melt_handler->add_current_constraints (new_current_constraints);
 
diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index 8a5a3742909..79a194d1b26 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -2271,17 +2271,61 @@ namespace aspect
 
   template <int dim>
   void
-  Simulator<dim>::replace_outflow_boundary_ids(const unsigned int offset)
+  Simulator<dim>::replace_outflow_boundary_ids(const unsigned int offset,
+                                               const bool is_composition,
+                                               const unsigned int composition_index)
   {
     const Quadrature<dim-1> &quadrature_formula = introspection.face_quadratures.temperature;
 
+    bool consider_darcy_velocity = false;
+    bool consider_fluid_velocity = false;
+    unsigned int porosity_idx = numbers::invalid_unsigned_int;
+
+    // Check to see if we are attempting to replace the boundary conditions for the temperature, or for the
+    // composition. If we are attempting to replace the composition boundary conditions, check to see if the current
+    // compositional field is advected with the darcy velocity (via the darcy field advection method) or with the
+    // fluid velocity (via the melt field advection method, or via a compositional field named porosity when melt
+    // transport is included).
+    if (is_composition)
+      {
+        if (introspection.compositional_field_methods[composition_index] == Parameters<dim>::AdvectionFieldMethod::fem_darcy_field)
+          {
+            consider_darcy_velocity = true;
+            porosity_idx = introspection.find_composition_type(CompositionalFieldDescription::porosity);
+          }
+        if (introspection.compositional_field_methods[composition_index] == Parameters<dim>::AdvectionFieldMethod::fem_melt_field ||
+            (parameters.include_melt_transport && introspection.get_indices_for_fields_of_type(CompositionalFieldDescription::porosity)[0] == composition_index) )
+          consider_fluid_velocity = true;
+      }
+
+    // If the compositional field uses the darcy velocity, we need to update the gradients, otherwise we do not
+    const UpdateFlags update_flags
+      = UpdateFlags(
+          consider_darcy_velocity
+          ?
+          update_values |
+          update_gradients |
+          update_normal_vectors |
+          update_quadrature_points |
+          update_JxW_values
+          :
+          update_values |
+          update_normal_vectors |
+          update_quadrature_points |
+          update_JxW_values);
+
     FEFaceValues<dim> fe_face_values (*mapping,
                                       finite_element,
                                       quadrature_formula,
-                                      update_values   | update_normal_vectors |
-                                      update_quadrature_points | update_JxW_values);
+                                      update_flags);
 
     std::vector<Tensor<1,dim>> face_current_velocity_values (fe_face_values.n_quadrature_points);
+    std::vector<Tensor<1,dim>> face_current_fluid_velocity_values (fe_face_values.n_quadrature_points);
+
+    MaterialModel::MaterialModelInputs<dim> in(fe_face_values.n_quadrature_points, introspection.n_compositional_fields);
+    MaterialModel::MaterialModelOutputs<dim> out(fe_face_values.n_quadrature_points, introspection.n_compositional_fields);
+    MeltHandler<dim>::create_material_model_outputs(out);
+    MaterialModel::MeltOutputs<dim> *fluid_out = out.template get_additional_output<MaterialModel::MeltOutputs<dim>>();
 
     const auto &tangential_velocity_boundaries =
       boundary_velocity_manager.get_tangential_boundary_velocity_indicators();
@@ -2311,11 +2355,23 @@ namespace aspect
                 fe_face_values.reinit (cell, face_number);
                 fe_face_values[introspection.extractors.velocities].get_function_values(current_linearization_point,
                                                                                         face_current_velocity_values);
-
                 // ... check if the face is an outflow boundary by integrating the normal velocities
                 // (flux through the boundary) as: int u*n ds = Sum_q u(x_q)*n(x_q) JxW(x_q)...
+                // do this for the solid velocity, the darcy velocity, or the fluid velocity.
                 double integrated_flow = 0;
 
+                if (consider_darcy_velocity)
+                  {
+                    in.reinit(fe_face_values, cell, introspection, solution);
+                    material_model->evaluate(in, out);
+                  }
+
+                if (consider_fluid_velocity)
+                  {
+                    const FEValuesExtractors::Vector ex_u_f = introspection.variable("fluid velocity").extractor_vector();
+                    fe_face_values[ex_u_f].get_function_values(current_linearization_point, face_current_fluid_velocity_values);
+                  }
+
                 for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
                   {
                     Tensor<1,dim> boundary_velocity;
@@ -2325,8 +2381,32 @@ namespace aspect
                       boundary_velocity = boundary_velocity_manager.boundary_velocity(face->boundary_id(),
                                                                                       fe_face_values.quadrature_point(q));
 
-                    integrated_flow += (boundary_velocity * fe_face_values.normal_vector(q)) *
-                                       fe_face_values.JxW(q);
+                    if (consider_darcy_velocity)
+                      {
+                        const double porosity = std::max(in.composition[q][porosity_idx], 1e-10);
+                        const Tensor<1,dim> gravity = gravity_model->gravity_vector(in.position[q]);
+                        const double solid_density = out.densities[q];
+                        const double fluid_viscosity = fluid_out->fluid_viscosities[q];
+                        const double fluid_density = fluid_out->fluid_densities[q];
+                        const double permeability = fluid_out->permeabilities[q];
+                        const Tensor<1,dim> boundary_darcy_velocity = boundary_velocity -
+                                                                      permeability / fluid_viscosity / porosity * gravity *
+                                                                      (solid_density - fluid_density);
+                        integrated_flow += (boundary_darcy_velocity * fe_face_values.normal_vector(q)) *
+                                           fe_face_values.JxW(q);
+                      }
+
+                    else if (consider_fluid_velocity)
+                      {
+                        integrated_flow += (face_current_fluid_velocity_values[q] * fe_face_values.normal_vector(q)) *
+                                           fe_face_values.JxW(q);
+                      }
+
+                    else
+                      {
+                        integrated_flow += (boundary_velocity * fe_face_values.normal_vector(q)) *
+                                           fe_face_values.JxW(q);
+                      }
                   }
 
                 // ... and change the boundary id of any outflow boundary faces.
@@ -2758,7 +2838,9 @@ namespace aspect
   template void Simulator<dim>::initialize_current_linearization_point (); \
   template void Simulator<dim>::interpolate_material_output_into_advection_field(const AdvectionField &adv_field); \
   template void Simulator<dim>::check_consistency_of_formulation(); \
-  template void Simulator<dim>::replace_outflow_boundary_ids(const unsigned int boundary_id_offset); \
+  template void Simulator<dim>::replace_outflow_boundary_ids(const unsigned int boundary_id_offset, \
+                                                             const bool is_composition, \
+                                                             const unsigned int composition_index); \
   template void Simulator<dim>::restore_outflow_boundary_ids(const unsigned int boundary_id_offset); \
   template void Simulator<dim>::check_consistency_of_boundary_conditions() const; \
   template double Simulator<dim>::compute_initial_newton_residual(const LinearAlgebra::BlockVector &linearized_stokes_initial_guess); \