diff --git a/source/material_model/reactive_fluid_transport.cc b/source/material_model/reactive_fluid_transport.cc
index c719593690c..13531e492b1 100644
--- a/source/material_model/reactive_fluid_transport.cc
+++ b/source/material_model/reactive_fluid_transport.cc
@@ -86,25 +86,44 @@ namespace aspect
                     std::vector<double> &melt_fractions,
                     const MaterialModel::MaterialModelOutputs<dim> *out) const
     {
-      const std::shared_ptr<const MeltOutputs<dim>> fluid_out = out->template get_additional_output_object<MeltOutputs<dim>>();
-      if (fluid_out != nullptr && in.requests_property(MaterialProperties::additional_outputs))
+      // the Katz 2003 model does not need information about fluid content
+      if (fluid_solid_reaction_scheme == katz2003)
         {
+          for (unsigned int q=0; q<in.n_evaluation_points(); ++q)
+            melt_fractions[q] = katz2003_model.melt_fraction(in.temperature[q],
+                                                             this->get_adiabatic_conditions().pressure(in.position[q]));
+        }
+      else
+        {
+          Assert(out != nullptr,
+                 ExcMessage("The material model 'ReactiveFluidTransport' requires the material model "
+                            "outputs in order to compute melt fractions, but none were provided."));
+          Assert(out->template has_additional_output_object<MeltOutputs<dim>>(),
+                 ExcMessage("The material model 'ReactiveFluidTransport' requires melt material "
+                            "outputs to compute the melt fractions, but none were provided."));
+
+          const std::shared_ptr<const MeltOutputs<dim>> fluid_out = out->template get_additional_output_object<MeltOutputs<dim>>();
+
+          const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
+          const unsigned int bound_fluid_idx = (fluid_solid_reaction_scheme == zero_solubility)
+                                               ?
+                                               this->introspection().compositional_index_for_name("bound_fluid")
+                                               :
+                                               numbers::invalid_unsigned_int;
+
           for (unsigned int q=0; q<in.n_evaluation_points(); ++q)
             {
-              const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
+              // This function outputs the mass fraction of equilibrium free water at each quadrature point. However,
+              // the free water is a volume fraction, while the bound water is a mass fraction. Convert the free water
+              // to a mass fraction so that we correctly compare the two values.
+              const double volume_fraction_porosity = in.composition[q][porosity_idx];
+              const double bulk_density = compute_bulk_density(volume_fraction_porosity, out->densities[q], fluid_out->fluid_densities[q]);
+              const double mass_fraction_porosity = compute_mass_fraction(volume_fraction_porosity, fluid_out->fluid_densities[q], bulk_density);
+
               switch (fluid_solid_reaction_scheme)
                 {
                   case no_reaction:
                   {
-                    Assert(out != nullptr,
-                           ExcMessage("The MaterialModelOutputs provided to the No Reactions reaction model is a nullptr."));
-                    const double volume_fraction_porosity = in.composition[q][porosity_idx];
-                    // Melt fractions outputs the mass fraction of equilibrium free water at each quadrature point. However,
-                    // the free water is a volume fraction, while the bound water is a mass fraction. Convert the free water
-                    // to a mass fraction so that we correctly compare the two values.
-                    const std::shared_ptr<const MeltOutputs<dim>> fluid_out = out->template get_additional_output_object<MeltOutputs<dim>>();
-                    const double bulk_density = compute_bulk_density(volume_fraction_porosity, out->densities[q], fluid_out->fluid_densities[q]);
-                    const double mass_fraction_porosity = compute_mass_fraction(volume_fraction_porosity, fluid_out->fluid_densities[q], bulk_density);
                     // No reactions occur between the solid and fluid phases,
                     // and the fluid mass fraction (stored in the melt_fractions
                     // vector) is equal to the mass fraction of the porosity.
@@ -116,39 +135,14 @@ namespace aspect
                     // The fluid volume fraction in equilibrium with the solid
                     // at any point (stored in the melt_fractions vector) is
                     // equal to the sum of the bound fluid content and porosity.
-                    Assert(out != nullptr,
-                           ExcMessage("The MaterialModelOutputs provided to the Zero Solubility reaction model is a nullptr."));
-                    const double volume_fraction_porosity = in.composition[q][porosity_idx];
-                    const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
-                    // Melt fractions outputs the mass fraction of equilibrium free water at each quadrature point. However,
-                    // the free water is a volume fraction, while the bound water is a mass fraction. Convert the free water
-                    // to a mass fraction so that we correctly compare the two values.
-                    const std::shared_ptr<const MeltOutputs<dim>> fluid_out = out->template get_additional_output_object<MeltOutputs<dim>>();
-                    const double bulk_density = compute_bulk_density(volume_fraction_porosity, out->densities[q], fluid_out->fluid_densities[q]);
-                    const double mass_fraction_porosity = compute_mass_fraction(volume_fraction_porosity, fluid_out->fluid_densities[q], bulk_density);
                     melt_fractions[q] = in.composition[q][bound_fluid_idx] + mass_fraction_porosity;
                     break;
                   }
                   case tian_approximation:
                   {
-                    Assert(out != nullptr,
-                           ExcMessage("The MaterialModelOutputs provided to the Tian 2019 reaction model is a nullptr."));
-                    // The tian2019_model melt_fraction function requires the mass fraction of free water at each quadrature
-                    // point as an input. However, the free water is stored in the porosity field as a volume fraction.
-                    // We therefore convert the free water to a mass fraction so that we have the correct input.
-                    const double volume_fraction_porosity = in.composition[q][porosity_idx];
-                    const std::shared_ptr<const MeltOutputs<dim>> fluid_out = out->template get_additional_output_object<MeltOutputs<dim>>();
-                    const double bulk_density = compute_bulk_density(volume_fraction_porosity, out->densities[q], fluid_out->fluid_densities[q]);;
-                    const double mass_fraction_porosity = compute_mass_fraction(volume_fraction_porosity, fluid_out->fluid_densities[q], bulk_density);
                     melt_fractions[q] = tian2019_model.melt_fraction(in, mass_fraction_porosity, q);
                     break;
                   }
-                  case katz2003:
-                  {
-                    melt_fractions[q] = katz2003_model.melt_fraction(in.temperature[q],
-                                                                     this->get_adiabatic_conditions().pressure(in.position[q]));
-                    break;
-                  }
                   default:
                   {
                     AssertThrow(false, ExcNotImplemented());
@@ -236,6 +230,10 @@ namespace aspect
           if (this->get_parameters().use_operator_splitting && reaction_rate_out != nullptr
               && in.requests_property(MaterialProperties::reaction_rates))
             {
+              Assert(fluid_out != nullptr,
+                     ExcMessage("The material model 'ReactiveFluidTransport' requires melt material "
+                                "outputs to compute reaction rates, but none were provided."));
+
               // Fill reaction rate outputs if the model uses operator splitting.
               // Specifically, change the porosity (representing the amount of free fluid)
               // based on the water solubility and the fluid content.