diff --git a/doc/modules/changes/20241201_lhy11009 b/doc/modules/changes/20241201_lhy11009
new file mode 100644
index 00000000000..86854255027
--- /dev/null
+++ b/doc/modules/changes/20241201_lhy11009
@@ -0,0 +1,6 @@
+Changed: Renamed the variable 'PhaseFunctionInputs::phase_index' to
+PhaseFunctionInputs::phase_transition_index'. The new variable name
+is more precise since it is used to index phase transitions rather than phases.
+Material models that make use of the old name will have to be adjusted.
+<br>
+(Haoyuan Li, 2024/12/01)
diff --git a/include/aspect/material_model/utilities.h b/include/aspect/material_model/utilities.h
index fb48807a162..e1b4c821cae 100644
--- a/include/aspect/material_model/utilities.h
+++ b/include/aspect/material_model/utilities.h
@@ -528,7 +528,7 @@ namespace aspect
                             const double pressure,
                             const double depth,
                             const double pressure_depth_derivative,
-                            const unsigned int phase_index);
+                            const unsigned int phase_transition_index);
 
         double temperature;
         double pressure;
@@ -546,7 +546,7 @@ namespace aspect
          * which only has information that there are two phase functions
          * and what their properties are.
          */
-        unsigned int phase_index;
+        unsigned int phase_transition_index;
       };
 
       /**
@@ -775,14 +775,14 @@ namespace aspect
 
           /**
            * Return the Clapeyron slope (dp/dT of the transition) for
-           * phase transition number @p phase_index.
+           * phase transition number @p phase_transition_index.
            */
-          double get_transition_slope (const unsigned int phase_index) const;
+          double get_transition_slope (const unsigned int phase_transition_index) const;
 
           /**
-           * Return the depth for phase transition number @p phase_index.
+           * Return the depth for phase transition number @p phase_transition_index.
            */
-          double get_transition_depth (const unsigned int phase_index) const;
+          double get_transition_depth (const unsigned int phase_transition_index) const;
 
           /**
            * Return how many phase transitions there are for each chemical composition.
diff --git a/source/material_model/grain_size.cc b/source/material_model/grain_size.cc
index 172b5625d12..1d4f9593cc6 100644
--- a/source/material_model/grain_size.cc
+++ b/source/material_model/grain_size.cc
@@ -576,7 +576,7 @@ namespace aspect
 
                   for (unsigned int k=0; k<n_phase_transitions; ++k)
                     {
-                      phase_inputs.phase_index = k;
+                      phase_inputs.phase_transition_index = k;
                       phase_function_values[k] = phase_function->compute_value(phase_inputs);
                     }
 
diff --git a/source/material_model/utilities.cc b/source/material_model/utilities.cc
index 0b8c1c3dfd0..271e71a0522 100644
--- a/source/material_model/utilities.cc
+++ b/source/material_model/utilities.cc
@@ -1149,14 +1149,14 @@ namespace aspect
                                                     const double pressure_,
                                                     const double depth_,
                                                     const double pressure_depth_derivative_,
-                                                    const unsigned int phase_index_)
+                                                    const unsigned int phase_transition_index_)
 
         :
         temperature(temperature_),
         pressure(pressure_),
         depth(depth_),
         pressure_depth_derivative(pressure_depth_derivative_),
-        phase_index(phase_index_)
+        phase_transition_index(phase_transition_index_)
       {}
 
 
@@ -1208,7 +1208,7 @@ namespace aspect
         unsigned int n_comp = 0;
         for (unsigned int n_relevant_fields  = 0 ; n_relevant_fields  < this->introspection().n_chemical_composition_fields() + 1 ; n_relevant_fields ++)
           {
-            if (in.phase_index < start_phase_transition_index + n_phase_transitions_per_chemical_composition[n_relevant_fields])
+            if (in.phase_transition_index < start_phase_transition_index + n_phase_transitions_per_chemical_composition[n_relevant_fields])
               {
                 n_comp = n_relevant_fields ;
                 break;
@@ -1246,7 +1246,7 @@ namespace aspect
           }
 
         // determine the value of phase function to facilitate the exact transition
-        if ((matched_phase_transition_index != numbers::invalid_unsigned_int) && in.phase_index <= matched_phase_transition_index)
+        if ((matched_phase_transition_index != numbers::invalid_unsigned_int) && in.phase_transition_index <= matched_phase_transition_index)
           function_value = 1.0;
         else
           function_value = 0.0;
@@ -1424,13 +1424,13 @@ namespace aspect
       double
       PhaseFunction<dim>::compute_value (const PhaseFunctionInputs<dim> &in) const
       {
-        AssertIndexRange (in.phase_index, transition_temperature_lower_limits.size());
-        AssertIndexRange (in.phase_index, transition_temperature_upper_limits.size());
+        AssertIndexRange (in.phase_transition_index, transition_temperature_lower_limits.size());
+        AssertIndexRange (in.phase_transition_index, transition_temperature_upper_limits.size());
 
         // the percentage of material that has undergone the transition
         double function_value;
-        if (in.temperature < transition_temperature_lower_limits[in.phase_index] ||
-            in.temperature >= transition_temperature_upper_limits[in.phase_index])
+        if (in.temperature < transition_temperature_lower_limits[in.phase_transition_index] ||
+            in.temperature >= transition_temperature_upper_limits[in.phase_transition_index])
           {
             // assign 0.0 if temperature is out of range
             function_value = 0.0;
@@ -1439,40 +1439,40 @@ namespace aspect
           {
             if (use_depth_instead_of_pressure)
               {
-                AssertIndexRange (in.phase_index, transition_depths.size());
+                AssertIndexRange (in.phase_transition_index, transition_depths.size());
 
                 // calculate the deviation from the transition point (convert temperature to depth)
-                double depth_deviation = in.depth - transition_depths[in.phase_index];
+                double depth_deviation = in.depth - transition_depths[in.phase_transition_index];
 
                 if (in.pressure_depth_derivative != 0.0)
                   {
-                    AssertIndexRange (in.phase_index, transition_slopes.size());
-                    AssertIndexRange (in.phase_index, transition_temperatures.size());
+                    AssertIndexRange (in.phase_transition_index, transition_slopes.size());
+                    AssertIndexRange (in.phase_transition_index, transition_temperatures.size());
 
-                    depth_deviation -= transition_slopes[in.phase_index] / in.pressure_depth_derivative
-                                       * (in.temperature - transition_temperatures[in.phase_index]);
+                    depth_deviation -= transition_slopes[in.phase_transition_index] / in.pressure_depth_derivative
+                                       * (in.temperature - transition_temperatures[in.phase_transition_index]);
                   }
 
                 // use delta function for width = 0
-                AssertIndexRange (in.phase_index, transition_widths.size());
-                if (transition_widths[in.phase_index] == 0)
+                AssertIndexRange (in.phase_transition_index, transition_widths.size());
+                if (transition_widths[in.phase_transition_index] == 0)
                   function_value = (depth_deviation > 0) ? 1. : 0.;
                 else
-                  function_value = 0.5*(1.0 + std::tanh(depth_deviation / transition_widths[in.phase_index]));
+                  function_value = 0.5*(1.0 + std::tanh(depth_deviation / transition_widths[in.phase_transition_index]));
               }
             else
               {
                 // calculate the deviation from the transition point (convert temperature to pressure)
-                AssertIndexRange (in.phase_index, transition_pressures.size());
-                const double pressure_deviation = in.pressure - transition_pressures[in.phase_index]
-                                                  - transition_slopes[in.phase_index] * (in.temperature - transition_temperatures[in.phase_index]);
+                AssertIndexRange (in.phase_transition_index, transition_pressures.size());
+                const double pressure_deviation = in.pressure - transition_pressures[in.phase_transition_index]
+                                                  - transition_slopes[in.phase_transition_index] * (in.temperature - transition_temperatures[in.phase_transition_index]);
 
                 // use delta function for width = 0
-                AssertIndexRange (in.phase_index, transition_pressure_widths.size());
-                if (transition_pressure_widths[in.phase_index] == 0)
+                AssertIndexRange (in.phase_transition_index, transition_pressure_widths.size());
+                if (transition_pressure_widths[in.phase_transition_index] == 0)
                   function_value = (pressure_deviation > 0) ? 1. : 0.;
                 else
-                  function_value = 0.5*(1.0 + std::tanh(pressure_deviation / transition_pressure_widths[in.phase_index]));
+                  function_value = 0.5*(1.0 + std::tanh(pressure_deviation / transition_pressure_widths[in.phase_transition_index]));
               }
           }
 
@@ -1499,30 +1499,30 @@ namespace aspect
         // phase transition based on depth
         if (use_depth_instead_of_pressure)
           {
-            const Point<dim,double> transition_point = this->get_geometry_model().representative_point(transition_depths[in.phase_index]);
-            const Point<dim,double> transition_plus_width = this->get_geometry_model().representative_point(transition_depths[in.phase_index] + transition_widths[in.phase_index]);
-            const Point<dim,double> transition_minus_width = this->get_geometry_model().representative_point(transition_depths[in.phase_index] - transition_widths[in.phase_index]);
+            const Point<dim,double> transition_point = this->get_geometry_model().representative_point(transition_depths[in.phase_transition_index]);
+            const Point<dim,double> transition_plus_width = this->get_geometry_model().representative_point(transition_depths[in.phase_transition_index] + transition_widths[in.phase_transition_index]);
+            const Point<dim,double> transition_minus_width = this->get_geometry_model().representative_point(transition_depths[in.phase_transition_index] - transition_widths[in.phase_transition_index]);
             transition_pressure = this->get_adiabatic_conditions().pressure(transition_point);
             pressure_width = 0.5 * (this->get_adiabatic_conditions().pressure(transition_plus_width)
                                     - this->get_adiabatic_conditions().pressure(transition_minus_width));
-            width_temp = transition_widths[in.phase_index];
+            width_temp = transition_widths[in.phase_transition_index];
           }
         // using pressure instead of depth to define the phase transition
         else
           {
-            transition_pressure = transition_pressures[in.phase_index];
-            pressure_width = transition_pressure_widths[in.phase_index];
-            width_temp = transition_pressure_widths[in.phase_index];
+            transition_pressure = transition_pressures[in.phase_transition_index];
+            pressure_width = transition_pressure_widths[in.phase_transition_index];
+            width_temp = transition_pressure_widths[in.phase_transition_index];
           }
 
         // calculate the deviation from the transition point
         const double pressure_deviation = in.pressure - transition_pressure
-                                          - transition_slopes[in.phase_index] * (in.temperature - transition_temperatures[in.phase_index]);
+                                          - transition_slopes[in.phase_transition_index] * (in.temperature - transition_temperatures[in.phase_transition_index]);
 
         // calculate the analytical derivative of the phase function
         if (
-          (in.temperature < transition_temperature_lower_limits[in.phase_index]) ||
-          (in.temperature >= transition_temperature_upper_limits[in.phase_index])
+          (in.temperature < transition_temperature_lower_limits[in.phase_transition_index]) ||
+          (in.temperature >= transition_temperature_upper_limits[in.phase_transition_index])
         )
           {
             // return 0 if temperature is out of range
@@ -1609,9 +1609,9 @@ namespace aspect
       template <int dim>
       double
       PhaseFunction<dim>::
-      get_transition_slope (const unsigned int phase_index) const
+      get_transition_slope (const unsigned int phase_transition_index) const
       {
-        return transition_slopes[phase_index];
+        return transition_slopes[phase_transition_index];
       }
 
 
@@ -1619,9 +1619,9 @@ namespace aspect
       template <int dim>
       double
       PhaseFunction<dim>::
-      get_transition_depth (const unsigned int phase_index) const
+      get_transition_depth (const unsigned int phase_transition_index) const
       {
-        return transition_depths[phase_index];
+        return transition_depths[phase_transition_index];
       }
 
 
diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc
index d76bac94313..47accfe3d27 100644
--- a/source/material_model/visco_plastic.cc
+++ b/source/material_model/visco_plastic.cc
@@ -80,7 +80,7 @@ namespace aspect
 
           for (unsigned int j=0; j < phase_function.n_phase_transitions(); ++j)
             {
-              phase_inputs.phase_index = j;
+              phase_inputs.phase_transition_index = j;
               phase_function_values[j] = phase_function.compute_value(phase_inputs);
             }
         }
@@ -141,7 +141,7 @@ namespace aspect
           // Compute value of phase functions
           for (unsigned int j=0; j < phase_function.n_phase_transitions(); ++j)
             {
-              phase_inputs.phase_index = j;
+              phase_inputs.phase_transition_index = j;
               phase_function_values[j] = phase_function.compute_value(phase_inputs);
             }
 
@@ -209,7 +209,7 @@ namespace aspect
                 {
                   for (unsigned int j=0; j < phase_function_discrete->n_phase_transitions(); ++j)
                     {
-                      phase_inputs.phase_index = j;
+                      phase_inputs.phase_transition_index = j;
                       phase_function_discrete_values[j] = phase_function_discrete->compute_value(phase_inputs);
                     }
                   isostrain_viscosities =
diff --git a/tests/composite_viscous_outputs_phases.cc b/tests/composite_viscous_outputs_phases.cc
index 3650747cfab..7415ac90e27 100644
--- a/tests/composite_viscous_outputs_phases.cc
+++ b/tests/composite_viscous_outputs_phases.cc
@@ -169,7 +169,7 @@ void f(const aspect::SimulatorAccess<dim> &simulator_access,
           // Compute value of phase functions
           for (unsigned int j=0; j < phase_function.n_phase_transitions(); ++j)
             {
-              phase_inputs.phase_index = j;
+              phase_inputs.phase_transition_index = j;
               phase_function_values[j] = phase_function.compute_value(phase_inputs);
             }