diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1f0cdbb1e1..04fb0edadf 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,8 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/).
 
 ### Added
 
+- MINOR Assertion on fields required for physical property models were added to the code. This is done to ensure that if a physical property model necessitating a certain field is employed and that field is not defined, the exception message thrown is comprehensible. [#1065](https://github.com/lethe-cfd/lethe/pull/1065)
+
 - MINOR A new `make_table_scalars_vectors` function has been added in `utilities.h`(`.cc`). [#1062](https://github.com/lethe-cfd/lethe/pull/1062)
 
 ## [Master] - 2024-03-11
diff --git a/include/core/density_model.h b/include/core/density_model.h
index bb48b6c97a..934efc88bf 100644
--- a/include/core/density_model.h
+++ b/include/core/density_model.h
@@ -94,18 +94,18 @@ class DensityConstant : public DensityModel
   /**
    * @brief Compute density value.
    *
-   * @param[in] fields_value Value of the various fields on which the property
-   * may depend.
+   * @param[in] field_values Value of the various fields on which the property
+   * may depend. The map stores a single value per field.
    *
-   * @note Here, the @p fields_value parameter is ignored since the density
+   * @note Here, the @p field_values parameter is ignored since the density
    * remains constant.
    *
    * @return Density value.
    */
   double
-  value(const std::map<field, double> &fields_value) override
+  value(const std::map<field, double> &field_values) override
   {
-    (void)fields_value;
+    (void)field_values;
     return density;
   }
 
@@ -113,7 +113,7 @@ class DensityConstant : public DensityModel
    * @brief Compute a vector of density values.
    *
    * @param[in] field_vectors Vectors of the fields on which the density may
-   * depend.
+   * depend. The map stores a vector of values per field.
    *
    * @param[out] property_vector Vector of computed density values.
    *
@@ -133,7 +133,7 @@ class DensityConstant : public DensityModel
    * respect to a specified field.
    *
    * @param[in] field_values Values of the various fields on which the density
-   * may depend.
+   * may depend. The map stores a single value per field.
    *
    * @param[in] id Indicator of the field with respect to which the jacobian
    * should be computed.
@@ -158,7 +158,7 @@ class DensityConstant : public DensityModel
    * fields.
    *
    * @param[in] field_vectors Vector of values of the fields used to evaluate
-   * the density.
+   * the density. The map stores a vector of values per field.
    *
    * @param[in] id Identifier of the field with respect to which a derivative
    * should be computed.
@@ -167,7 +167,7 @@ class DensityConstant : public DensityModel
    * with respect to the field of the specified @p id. In this case, it returns
    * a vector of zeros since the density remains constant.
    *
-   * @note Here, the @p field_values and @p id parameters are ignored since the
+   * @note Here, the @p field_vectors and @p id parameters are ignored since the
    * density remains constant.
    *
    */
@@ -256,15 +256,19 @@ class DensityIsothermalIdealGas : public DensityModel
   /**
    * @brief Compute the density of the isothermal ideal gas.
    *
-   * @param[in] fields_value Value of the various fields on which the property
-   * may depend. In this case, the density depends on pressure.
+   * @param[in] field_values Values of the various fields on which the property
+   * may depend. In this case, the density depends on pressure. The map stores a
+   * single value per field.
    *
-   * @return Value of the density computed with the @p fields_value.
+   * @return Value of the density computed with the @p field_values.
    */
   double
-  value(const std::map<field, double> &fields_value) override
+  value(const std::map<field, double> &field_values) override
   {
-    const double pressure = fields_value.at(field::pressure);
+    AssertThrow(field_values.find(field::pressure) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "DensityIsothermalIdealGas", "pressure"));
+    const double pressure = field_values.at(field::pressure);
     return density_ref + psi * pressure;
   }
 
@@ -272,7 +276,8 @@ class DensityIsothermalIdealGas : public DensityModel
    * @brief Compute a vector of density values for an isothermal ideal gas.
    *
    * @param[in] field_vectors Vectors of the fields on which the density may
-   * depend. In this case, the density depends on pressure.
+   * depend. In this case, the density depends on pressure. The map stores a
+   * vector of values per field.
    *
    * @param[out] property_vector Vectors of computed density values.
    */
@@ -280,6 +285,9 @@ class DensityIsothermalIdealGas : public DensityModel
   vector_value(const std::map<field, std::vector<double>> &field_vectors,
                std::vector<double> &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::pressure) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "DensityIsothermalIdealGas", "pressure"));
     const std::vector<double> &pressure = field_vectors.at(field::pressure);
     for (unsigned int i = 0; i < property_vector.size(); ++i)
       property_vector[i] = density_ref + psi * pressure[i];
@@ -290,7 +298,7 @@ class DensityIsothermalIdealGas : public DensityModel
    * respect to a specified field.
    *
    * @param[in] field_values Values of the various fields on which the density
-   * may depend.
+   * may depend. The map stores a single value per field.
    *
    * @param[in] id Indicator of the field with respect to which the jacobian
    * should be computed.
@@ -313,7 +321,7 @@ class DensityIsothermalIdealGas : public DensityModel
    * an isothermal ideal gas.
    *
    * @param[in] field_vectors Vector of values of the fields used to evaluate
-   * the density.
+   * the density. The map stores a vector of values per field.
    *
    * @param[in] id Identifier of the field with respect to which a derivative
    * should be computed.
diff --git a/include/core/evaporation_model.h b/include/core/evaporation_model.h
index 1c98bc7e67..322025a09e 100644
--- a/include/core/evaporation_model.h
+++ b/include/core/evaporation_model.h
@@ -215,7 +215,7 @@ class EvaporationModelConstant : public EvaporationModel
     , ambient_gas_density_inv(1.0 / p_evaporation.ambient_gas_density)
     , liquid_density_inv(1.0 / p_evaporation.liquid_density)
   {
-    model_depends_on[temperature] = false;
+    model_depends_on[field::temperature] = false;
   }
 
 
@@ -413,7 +413,7 @@ class EvaporationModelTemperature : public EvaporationModel
     , liquid_density_inv(1.0 / p_evaporation.liquid_density)
     , universal_gas_constant(p_evaporation.universal_gas_constant)
   {
-    model_depends_on[temperature] = true;
+    model_depends_on[field::temperature] = true;
   }
 
   /**
@@ -425,6 +425,9 @@ class EvaporationModelTemperature : public EvaporationModel
   double
   saturation_pressure(const std::map<field, double> &field_values)
   {
+    AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "EvaporationModelTemperature", "temperature"));
     const double temperature_inv =
       1.0 / (field_values.at(field::temperature) + 1e-16);
 
@@ -450,6 +453,9 @@ class EvaporationModelTemperature : public EvaporationModel
   saturation_pressure(const std::map<field, std::vector<double>> &field_vectors,
                       std::vector<double> &saturation_pressure_vector)
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "EvaporationModelTemperature", "temperature"));
     const std::vector<double> &temperature =
       field_vectors.at(field::temperature);
 
@@ -506,6 +512,9 @@ class EvaporationModelTemperature : public EvaporationModel
   mass_flux(const std::map<field, std::vector<double>> &field_vectors,
             std::vector<double> &mass_flux_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "EvaporationModelTemperature", "temperature"));
     const std::vector<double> &temperature =
       field_vectors.at(field::temperature);
 
@@ -578,6 +587,9 @@ class EvaporationModelTemperature : public EvaporationModel
   heat_flux_jacobian(const std::map<field, double> &field_values,
                      const field                    id) override
   {
+    AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "EvaporationModelTemperature", "temperature"));
     const double temperature_inv =
       1.0 / (field_values.at(field::temperature) + 1e-16);
 
@@ -620,6 +632,10 @@ class EvaporationModelTemperature : public EvaporationModel
 
     if (id == field::temperature)
       {
+        AssertThrow(
+          field_vectors.find(field::temperature) != field_vectors.end(),
+          PhysicialPropertyModelFieldUndefined("EvaporationModelTemperature",
+                                               "temperature"));
         const std::vector<double> &temperature =
           field_vectors.at(field::temperature);
 
diff --git a/include/core/mobility_cahn_hilliard_model.h b/include/core/mobility_cahn_hilliard_model.h
index a72342db1e..91a2dbcaeb 100644
--- a/include/core/mobility_cahn_hilliard_model.h
+++ b/include/core/mobility_cahn_hilliard_model.h
@@ -199,6 +199,10 @@ class MobilityCahnHilliardModelQuartic : public MobilityCahnHilliardModel
   double
   value(const std::map<field, double> &fields_value) override
   {
+    AssertThrow(
+      fields_value.find(field::phase_order_cahn_hilliard) != fields_value.end(),
+      PhysicialPropertyModelFieldUndefined("MobilityCahnHilliardModelQuartic",
+                                           "phase_order_cahn_hilliard"));
     const double &phase_order_cahn_hilliard =
       fields_value.at(field::phase_order_cahn_hilliard);
 
@@ -221,6 +225,11 @@ class MobilityCahnHilliardModelQuartic : public MobilityCahnHilliardModel
   vector_value(const std::map<field, std::vector<double>> &field_vectors,
                std::vector<double> &property_vector) override
   {
+    AssertThrow(
+      field_vectors.find(field::phase_order_cahn_hilliard) !=
+        field_vectors.end(),
+      PhysicialPropertyModelFieldUndefined("MobilityCahnHilliardModelQuartic",
+                                           "phase_order_cahn_hilliard"));
     const std::vector<double> &phase_order_cahn_hilliard =
       field_vectors.at(field::phase_order_cahn_hilliard);
     for (unsigned int i = 0; i < property_vector.size(); ++i)
@@ -248,6 +257,10 @@ class MobilityCahnHilliardModelQuartic : public MobilityCahnHilliardModel
   double
   jacobian(const std::map<field, double> &fields_value, field /*id*/) override
   {
+    AssertThrow(
+      fields_value.find(field::phase_order_cahn_hilliard) != fields_value.end(),
+      PhysicialPropertyModelFieldUndefined("MobilityCahnHilliardModelQuartic",
+                                           "phase_order_cahn_hilliard"));
     const double &phase_order_cahn_hilliard =
       fields_value.at(field::phase_order_cahn_hilliard);
 
@@ -275,6 +288,11 @@ class MobilityCahnHilliardModelQuartic : public MobilityCahnHilliardModel
                   const field /*id*/,
                   std::vector<double> &jacobian_vector) override
   {
+    AssertThrow(
+      field_vectors.find(field::phase_order_cahn_hilliard) !=
+        field_vectors.end(),
+      PhysicialPropertyModelFieldUndefined("MobilityCahnHilliardModelQuartic",
+                                           "phase_order_cahn_hilliard"));
     const std::vector<double> &phase_order_cahn_hilliard =
       field_vectors.at(field::phase_order_cahn_hilliard);
     for (unsigned int i = 0; i < jacobian_vector.size(); ++i)
diff --git a/include/core/parameters.h b/include/core/parameters.h
index b505732d89..d90ee2d777 100644
--- a/include/core/parameters.h
+++ b/include/core/parameters.h
@@ -254,7 +254,7 @@ namespace Parameters
   };
 
   /**
-   * @brief Carreau rheological model to solve for non Newtonian
+   * @brief Carreau model to solve for non Newtonian
    * flows.
    */
   struct CarreauParameters
diff --git a/include/core/physical_property_model.h b/include/core/physical_property_model.h
index ce8ae53633..f6dbf02720 100644
--- a/include/core/physical_property_model.h
+++ b/include/core/physical_property_model.h
@@ -33,6 +33,13 @@ DeclException2(SizeOfFields,
                << " is not equal to the number of values for another field "
                << arg2);
 
+DeclException2(PhysicialPropertyModelFieldUndefined,
+               std::string,
+               std::string,
+               << "Error in '" << arg1 << "' model. \n "
+               << "The '" << arg2
+               << "' field required by the model is not defined.");
+
 /*
  * Fields on which physical property can depend. All fields are assumed
  * to be at time t+dt other than those for which a _p suffix is explicitly
diff --git a/include/core/rheological_model.h b/include/core/rheological_model.h
index 37f62e4999..b848fca6e4 100644
--- a/include/core/rheological_model.h
+++ b/include/core/rheological_model.h
@@ -300,8 +300,8 @@ class PowerLaw : public RheologicalModel
     , n(n)
     , shear_rate_min(shear_rate_min)
   {
-    this->model_depends_on[shear_rate]    = true;
-    this->non_newtonian_rheological_model = true;
+    this->model_depends_on[field::shear_rate] = true;
+    this->non_newtonian_rheological_model     = true;
   }
 
   /**
@@ -404,6 +404,8 @@ class PowerLaw : public RheologicalModel
     const double                  &p_density_ref,
     const std::map<field, double> &field_values) const override
   {
+    AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined("PowerLaw", "shear_rate"));
     const double shear_rate_magnitude = field_values.at(field::shear_rate);
     return calculate_kinematic_viscosity(shear_rate_magnitude) * p_density_ref;
   }
@@ -420,6 +422,8 @@ class PowerLaw : public RheologicalModel
     const std::map<field, std::vector<double>> &field_vectors,
     std::vector<double>                        &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined("PowerLaw", "shear_rate"));
     const std::vector<double> &shear_rate_magnitude =
       field_vectors.at(field::shear_rate);
 
@@ -469,8 +473,8 @@ class Carreau : public RheologicalModel
     , a(p_a)
     , n(p_n)
   {
-    this->model_depends_on[shear_rate]    = true;
-    this->non_newtonian_rheological_model = true;
+    this->model_depends_on[field::shear_rate] = true;
+    this->non_newtonian_rheological_model     = true;
   }
 
   /**
@@ -564,6 +568,8 @@ class Carreau : public RheologicalModel
     const double                  &p_density_ref,
     const std::map<field, double> &field_values) const override
   {
+    AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined("Carreau", "shear_rate"));
     const double shear_rate_magnitude = field_values.at(field::shear_rate);
     return calculate_kinematic_viscosity(shear_rate_magnitude) * p_density_ref;
   }
@@ -580,6 +586,8 @@ class Carreau : public RheologicalModel
     const std::map<field, std::vector<double>> &field_vectors,
     std::vector<double>                        &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined("Carreau", "shear_rate"));
     const std::vector<double> &shear_rate_magnitude =
       field_vectors.at(field::shear_rate);
 
@@ -639,7 +647,7 @@ class PhaseChangeRheology : public RheologicalModel
   PhaseChangeRheology(const Parameters::PhaseChange p_phase_change_parameters)
     : param(p_phase_change_parameters)
   {
-    this->model_depends_on[temperature] = true;
+    this->model_depends_on[field::temperature] = true;
   }
 
   /**
@@ -700,6 +708,9 @@ class PhaseChangeRheology : public RheologicalModel
     const double                  &p_density_ref,
     const std::map<field, double> &field_values) const override
   {
+    AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                     "temperature"));
     const double temperature = field_values.at(field::temperature);
     return kinematic_viscosity(temperature) * p_density_ref;
   }
@@ -716,6 +727,9 @@ class PhaseChangeRheology : public RheologicalModel
     const std::map<field, std::vector<double>> &field_vectors,
     std::vector<double>                        &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                     "temperature"));
     const std::vector<double> &temperature =
       field_vectors.at(field::temperature);
 
diff --git a/include/core/specific_heat_model.h b/include/core/specific_heat_model.h
index 3054e91795..513ee03383 100644
--- a/include/core/specific_heat_model.h
+++ b/include/core/specific_heat_model.h
@@ -25,20 +25,20 @@
 using namespace dealii;
 
 /**
- * @brief Abstract class that allows to calculate the
- * specific heat on each quadrature point using the temperature of the fluid.
- * magnitude. SpecficiHeat::get_specific_heat() is a pure virtual method,
- * since it can only be calculated knowing the model for the specific
- * heat that has been specifid
+ * @brief Abstract class for calculating the specific heat of materials on
+ * quadrature points using the temperature of the fluid.
  */
 class SpecificHeatModel : public PhysicalPropertyModel
 {
 public:
   /**
-   * @brief Instantiates and returns a pointer to a SpecificHeatModel object by casting it to
-   * the proper child class
+   * @brief Instantiates and returns a pointer to a SpecificHeatModel object
+   * by casting it to the proper child class.
    *
-   * @param material_properties Parameters for a material
+   * @param[in] material_properties Property parameters of a material (fluid or
+   * solid).
+   *
+   * @return Casted SpecificHeatModel object.
    */
   static std::shared_ptr<SpecificHeatModel>
   model_cast(const Parameters::Material &material_properties);
@@ -46,68 +46,108 @@ class SpecificHeatModel : public PhysicalPropertyModel
 
 
 /**
- * @brief Constant specific heat. Returns a constant specific
- * heat for a fluid
+ * @brief Constant specific heat model.
  */
 class ConstantSpecificHeat : public SpecificHeatModel
 {
 public:
   /**
-   * @brief Default constructor
+   * @brief Constructor of the constant specific heat model.
+   *
+   * @param[in] p_specific_heat Constant specific heat value.
    */
   ConstantSpecificHeat(const double p_specific_heat)
     : specific_heat(p_specific_heat)
   {}
 
   /**
-   * @brief value Calculates the value of the specific heat.
-   * @param fields_value Value of the various field on which the specific heat depends.
-   * @return value of the specific heat.
+   * @brief Compute specific heat value.
+   *
+   * @param[in] field_values Value of the various field on which the specific
+   * heat may depend. The map stores a single value per field.
+   *
+   * @note Here, the @p field_values parameter is ignored since the specific
+   * heat remains constant.
+   *
+   * @return Specific heat value.
    */
   double
-  value(const std::map<field, double> & /*fields_value*/) override
+  value(const std::map<field, double> &field_values) override
   {
+    (void)field_values;
     return specific_heat;
   };
 
   /**
-   * @brief vector_value Calculates the vector of specific heat.
-   * @param field_vectors Vector of fields on which the specific heat may depend.
-   * @param property_vector Vector of specific_heat values.
+   * @brief vector_value Compute specific heat values.
+   *
+   * @param[in] field_vectors Vectors of the fields on which the specific heat
+   * may depend. The map stores a vector of values per field.
+   *
+   * @param[out] property_vector Vector of computed specific heat values.
+   *
+   * @note Here, the @p field_vectors parameter is ignored since the specific
+   * heat remains constant.
    */
   void
-  vector_value(const std::map<field, std::vector<double>> & /*field_vectors*/,
+  vector_value(const std::map<field, std::vector<double>> &field_vectors,
                std::vector<double> &property_vector) override
   {
+    (void)field_vectors;
     std::fill(property_vector.begin(), property_vector.end(), specific_heat);
   }
 
   /**
-   * @brief jacobian Calculates the jacobian (the partial derivative) of the specific heat with respect to a field
-   * @param field_values Value of the various fields on which the specific heat may depend.
-   * @param id Indicator of the field with respect to which the jacobian
-   * should be calculated
-   * @return value of the partial derivative of the specific heat with respect to the field.
+   * @brief Compute the jacobian (the partial derivative) of the specific heat
+   * with respect to a specified field.
+   *
+   * @param[in] field_values Values of the various fields on which the specific
+   * heat may depend. The map stores a single value per field.
+   *
+   * @param[in] id Indicator of the field with respect to which the jacobian
+   * should be computed.
+   *
+   * @note Here, the @p field_values and @p id parameters are ignored since the
+   * specific heat remains constant.
+   *
+   * @return Value of the derivative of the specific heat with respect to the
+   * specified field. Since the specific heat remains constant, this function
+   * returns zero.
    */
   double
-  jacobian(const std::map<field, double> & /*field_values*/,
-           field /*id*/) override
+  jacobian(const std::map<field, double> &field_values, field id) override
   {
+    (void)field_values;
+    (void)id;
     return 0;
   };
 
   /**
-   * @brief vector_jacobian Calculate the derivative of the specific heat with respect to a field
-   * @param field_vectors Vector for the values of the fields used to evaluate the property
-   * @param id Identifier of the field with respect to which a derivative should be calculated
-   * @param jacobian Vector of the value of the derivative of the specific heat with respect to the field id
+   * @brief Computes the derivative of the specific heat with respect to
+   * specified fields.
+   *
+   * @param[in] field_vectors Vector of values of the fields used to evaluate
+   * the specific heat. The map stores a vector of values per field.
+   *
+   * @param[in] id Identifier of the field with respect to which a derivative
+   * should be computed.
+   *
+   * @param[out] jacobian Vector of computed derivative values of the specific
+   * heat
+   * with respect to the field of the specified @p id. In this case, it returns
+   * a vector of zeros since the specific heat remains constant.
+   *
+   * @note Here, the @p field_vectors and @p id parameters are ignored since the
+   * specific heat remains constant.
+   *
    */
   void
-  vector_jacobian(
-    const std::map<field, std::vector<double>> & /*field_vectors*/,
-    const field /*id*/,
-    std::vector<double> &jacobian_vector) override
+  vector_jacobian(const std::map<field, std::vector<double>> &field_vectors,
+                  const field                                 id,
+                  std::vector<double> &jacobian_vector) override
   {
+    (void)field_vectors;
+    (void)id;
     std::fill(jacobian_vector.begin(), jacobian_vector.end(), 0);
   };
 
@@ -117,15 +157,35 @@ class ConstantSpecificHeat : public SpecificHeatModel
 
 
 /**
- * @brief This model takes into account the phase change of a material
- * by considering the latent heat into the specific heat over a phase change
- * interval determined by [T_solidus,T_liquidus].
+ * @brief Phase change specific heat model.
+ *
+ * The phase change specific heat \f$(c^{*}_\mathrm{p})\f$ model takes into
+ * account the phase change of a material by considering the latent heat into
+ * the specific heat over a phase change temperature interval determined by
+ * \f$\left[T_\mathrm{s},T_\mathrm{l}\right]\f$ where \f$T_\mathrm{s}\f$ is the
+ * solidus temperature and \f$T_\mathrm{l}\f$ is the liquidus temperature.
+ *
+ * The specific heat is evaluated as:
+ * \f$ c^{*}_\mathrm{p}(T)  =
+ *    \begin{cases}
+ *      C_\mathrm{p,s} & \text{if} \; T<T_\mathrm{s}\\
+ *      \frac{C_\mathrm{p,s}+C_\mathrm{p,l}}{2}+\frac{h_\mathrm{l}}
+ *        {T_\mathrm{l}-T_\mathrm{s}} & \text{if} \;
+ * T\in[T_\mathrm{s},T_\mathrm{l}]\\ C_\mathrm{p,l} & \text{if} \;
+ * T>T_\mathrm{l} \end{cases} \f$
+ *
+ * where \f$C_\mathrm{p,s}\f$ and \f$C_\mathrm{p,l}\f$ are the solid and liquid
+ * specific heat, respectively, and \f$h_\mathrm{l}\f$ is the latent enthalpy
+ * (enthalpy related to the phase change).
  */
 class PhaseChangeSpecificHeat : public SpecificHeatModel
 {
 public:
   /**
-   * @brief Default constructor
+   * @brief Constructor of the phase change specific heat model.
+   *
+   * @param[in] p_phase_change_parameters Set of parameters of the phase change
+   * model.
    */
   PhaseChangeSpecificHeat(
     const Parameters::PhaseChange p_phase_change_parameters)
@@ -137,15 +197,29 @@ class PhaseChangeSpecificHeat : public SpecificHeatModel
   }
 
   /**
-   * @brief value Calculates the value of the phase change specific heat.
-   * @param fields_value Value of the various field on which the specific heat depends.
-   * @return value of the specific heat.   */
+   * @brief Compute the specific heat of the isothermal ideal gas.
+   *
+   * @param[in] field_values Values of the various fields on which the property
+   * may depend. In this case, the specific heat depends on temperature.
+   * The map stores a single value per field.
+   *
+   * @return Value of the specific heat computed with the @p field_values.
+   */
   double
-  value(const std::map<field, double> &fields_value) override
+  value(const std::map<field, double> &field_values) override
   {
-    double temperature    = fields_value.at(field::temperature);
-    double temperature_p1 = fields_value.at(field::temperature_p1);
-    double temperature_p2 = fields_value.at(field::temperature_p2);
+    AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined("PhaseChangeSpecificHeat",
+                                                     "temperature"));
+    AssertThrow(field_values.find(field::temperature_p1) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined("PhaseChangeSpecificHeat",
+                                                     "temperature_p1"));
+    AssertThrow(field_values.find(field::temperature_p2) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined("PhaseChangeSpecificHeat",
+                                                     "temperature_p2"));
+    double temperature    = field_values.at(field::temperature);
+    double temperature_p1 = field_values.at(field::temperature_p1);
+    double temperature_p2 = field_values.at(field::temperature_p2);
 
 
     // Gather information required from the simulation control to have the time
@@ -176,8 +250,6 @@ class PhaseChangeSpecificHeat : public SpecificHeatModel
             const double dT =
               bdf_coefs[0] * temperature + bdf_coefs[1] * temperature_p1;
             return dH / dT;
-
-            break;
           }
 
         case Parameters::SimulationControl::TimeSteppingMethod::bdf2:
@@ -198,7 +270,6 @@ class PhaseChangeSpecificHeat : public SpecificHeatModel
                               bdf_coefs[1] * temperature_p1 +
                               bdf_coefs[2] * temperature_p2;
             return dH / dT;
-            break;
           }
 
         default:
@@ -207,16 +278,30 @@ class PhaseChangeSpecificHeat : public SpecificHeatModel
       }
   }
 
-
   /**
-   * @brief vector_value Calculates the vector of specific heat.
-   * @param field_vectors Vector of fields on which the specific heat may depend.
-   * @param property_vector Vector of specific_heat values.
+   * @brief Compute a vector of specific heat values for an isothermal ideal gas.
+   *
+   * @param[in] field_vectors Vectors of the fields on which the specific heat
+   * may depend. In this case, the specific heat depends on temperature. The map
+   * stores a vector of values per field.
+   *
+   * @param[out] property_vector Vectors of computed specific heat values.
    */
   void
   vector_value(const std::map<field, std::vector<double>> &field_vectors,
                std::vector<double> &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined("PhaseChangeSpecificHeat",
+                                                     "temperature"));
+    AssertThrow(field_vectors.find(field::temperature_p1) !=
+                  field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined("PhaseChangeSpecificHeat",
+                                                     "temperature_p1"));
+    AssertThrow(field_vectors.find(field::temperature_p2) !=
+                  field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined("PhaseChangeSpecificHeat",
+                                                     "temperature_p2"));
     const std::vector<double> &temperature_vec =
       field_vectors.at(field::temperature);
     const std::vector<double> &p1_temperature_vec =
@@ -300,55 +385,80 @@ class PhaseChangeSpecificHeat : public SpecificHeatModel
   }
 
   /**
-   * @brief jacobian Calculates the jacobian (the partial derivative) of the specific heat with respect to a field
-   * @param field_values Value of the various fields on which the specific heat may depend.
-   * @param id Indicator of the field with respect to which the jacobian
-   * should be calculated
-   * @return value of the partial derivative of the specific heat with respect to the field.
+   * @brief Compute the jacobian (the partial derivative) of the specific heat
+   * with respect to a specified field.
+   *
+   * @param[in] field_values Values of the various fields on which the specific
+   * heat may depend. The map stores a single value per field.
+   *
+   * @param[in] id Indicator of the field with respect to which the jacobian
+   * should be computed.
+   *
+   * @return Value of the derivative of the specific heat with respect to the
+   * specified field.
    */
   double
   jacobian(const std::map<field, double> &field_values, field id) override
   {
     if (id == field::temperature)
-      return numerical_jacobian(field_values, field::temperature);
+      {
+        AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                    PhysicialPropertyModelFieldUndefined(
+                      "EvaporationModelTemperature", "temperature"));
+        return numerical_jacobian(field_values, field::temperature);
+      }
     else
       return 0;
   };
 
   /**
-   * @brief vector_jacobian Calculate the derivative of the specific heat with respect to a field
-   * @param field_vectors Vector for the values of the fields used to evaluate the property
-   * @param id Identifier of the field with respect to which a derivative should be calculated
-   * @param jacobian Vector of the value of the derivative of the specific heat with respect to the field id
+   * @brief Compute the derivative of the specific heat with respect to a field
+   * for an isothermal ideal gas.
+   *
+   * @param[in] field_vectors Vector of values of the fields used to evaluate
+   * the specific heat. The map stores a vector of values per field.
+   *
+   * @param[in] id Identifier of the field with respect to which a derivative
+   * should be computed.
+   *
+   * @param[out] jacobian Vector of computed derivative values of the specific
+   * heat with respect to the field of the specified @p id.
    */
   void
   vector_jacobian(const std::map<field, std::vector<double>> &field_vectors,
                   const field                                 id,
                   std::vector<double> &jacobian_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "EvaporationModelTemperature", "temperature"));
     vector_numerical_jacobian(field_vectors, id, jacobian_vector);
   };
 
   /**
-   * @brief enthalpy Calculates the enthalpy of a phase change material for a temperature T
-   * The enthalpy is defined as :
+   * @brief Compute the enthalpy of a phase change material for a temperature
+   * \f$T\f$.
+   *
+   * The enthalpy \f$(H)\f$ is defined:
+   *
+   * - In pure liquid \f$(T>T_\mathrm{l})\f$ as:\n
+   * \f$H(T) = c_\mathrm{p,s} T_\mathrm{s} + 0.5(c_\mathrm{p,s}+c_\mathrm{p,l})
+   * (T_\mathrm{l}-T_\mathrm{s}) + h_\mathrm{l} + c_\mathrm{p,l}
+   * (T-T_\mathrm{l})\f$\n\n
    *
-   * Pure liquid
-   * ------------
-   * if (T>T_liquidus) : H = cp_solid * T_solidus + 0.5*(cp_solid+cp_liquid) *
-   * (T_liquidus-T_solidus) + latent_enthalpy + cp_liquid * (T-T_liquidus)
+   * - In liquid-solid mixture \f$(T_\mathrm{s} \leq T \leq T_\mathrm{l})\f$
+   * as:\n
+   *\f$H(T) = c_\mathrm{p,s} T_\mathrm{s} + 0.5(c_\mathrm{p,s}+c_\mathrm{p,l})
+   * (T_\mathrm{l}-T_\mathrm{s}) + \alpha_\mathrm{l} h_\mathrm{l}\f$\n
+   * where \f$\alpha_\mathrm{l} = \frac{T-T_\mathrm{s}}{T_\mathrm{l}-
+   * T_\mathrm{s}}\f$ is the liquid fraction.\n\n
    *
-   * Liquid-solid mixture
-   * -----------------
-   * else if (T>T_solidus) : cp_solid * T_solidus + 0.5*(cp_solid+cp_liquid) *
-   * (T-T_solidus) + liquid_fraction * latent_enthalpy
+   * - In pure solid \f$(T<T_\mathrm{s})\f$ as:\n
+   * \f$H(T) = c_\mathrm{p,s} T\f$
    *
-   * Pure solid
-   * ------------
-   * else : cp_solid * T
+   * @param[in] T temperature at which the enthalpy is computed.
    *
-   * @param T temperature at which to calculate the enthalpy
-   * @return Value of the enthalpy
+   * @return Value of the enthalpy evaluated at \f$T\f$.
    */
 
   inline double
diff --git a/include/core/surface_tension_model.h b/include/core/surface_tension_model.h
index 8d4b7fb306..54d9ca6f15 100644
--- a/include/core/surface_tension_model.h
+++ b/include/core/surface_tension_model.h
@@ -170,6 +170,9 @@ class SurfaceTensionLinear : public SurfaceTensionModel
   double
   value(const std::map<field, double> &fields_value) override
   {
+    AssertThrow(fields_value.find(field::temperature) != fields_value.end(),
+                PhysicialPropertyModelFieldUndefined("SurfaceTensionLinear",
+                                                     "temperature"));
     const double temperature = fields_value.at(field::temperature);
     return surface_tension_coefficient +
            surface_tension_gradient * (temperature - T_0);
@@ -185,6 +188,9 @@ class SurfaceTensionLinear : public SurfaceTensionModel
   vector_value(const std::map<field, std::vector<double>> &field_vectors,
                std::vector<double> &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined("SurfaceTensionLinear",
+                                                     "temperature"));
     const std::vector<double> &temperature =
       field_vectors.at(field::temperature);
     for (unsigned int i = 0; i < property_vector.size(); ++i)
@@ -273,6 +279,9 @@ class SurfaceTensionPhaseChange : public SurfaceTensionLinear
   double
   value(const std::map<field, double> &fields_value) override
   {
+    AssertThrow(fields_value.find(field::temperature) != fields_value.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "SurfaceTensionPhaseChange", "temperature"));
     const double temperature = fields_value.at(field::temperature);
 
     double surface_tension;
@@ -304,6 +313,9 @@ class SurfaceTensionPhaseChange : public SurfaceTensionLinear
   vector_value(const std::map<field, std::vector<double>> &field_vectors,
                std::vector<double> &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "SurfaceTensionPhaseChange", "temperature"));
     const std::vector<double> &temperature =
       field_vectors.at(field::temperature);
     for (unsigned int i = 0; i < property_vector.size(); ++i)
@@ -339,6 +351,9 @@ class SurfaceTensionPhaseChange : public SurfaceTensionLinear
   {
     if (id == field::temperature)
       {
+        AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                    PhysicialPropertyModelFieldUndefined(
+                      "SurfaceTensionPhaseChange", "temperature"));
         const double temperature = field_values.at(field::temperature);
         if (temperature < T_solidus)
           return 0;
@@ -372,6 +387,10 @@ class SurfaceTensionPhaseChange : public SurfaceTensionLinear
   {
     if (id == field::temperature)
       {
+        AssertThrow(
+          field_vectors.find(field::temperature) != field_vectors.end(),
+          PhysicialPropertyModelFieldUndefined("SurfaceTensionPhaseChange",
+                                               "temperature"));
         const std::vector<double> &temperature =
           field_vectors.at(field::temperature);
         for (unsigned int i = 0; i < jacobian_vector.size(); ++i)
diff --git a/include/core/thermal_conductivity_model.h b/include/core/thermal_conductivity_model.h
index 54c1b43fea..f9a16f561e 100644
--- a/include/core/thermal_conductivity_model.h
+++ b/include/core/thermal_conductivity_model.h
@@ -136,6 +136,9 @@ class ThermalConductivityLinear : public ThermalConductivityModel
   double
   value(const std::map<field, double> &fields_value) override
   {
+    AssertThrow(fields_value.find(field::temperature) != fields_value.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "ThermalConductivityLinear", "temperature"));
     return A + B * fields_value.at(field::temperature);
   };
 
@@ -148,6 +151,9 @@ class ThermalConductivityLinear : public ThermalConductivityModel
   vector_value(const std::map<field, std::vector<double>> &field_vectors,
                std::vector<double> &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "ThermalConductivityLinear", "temperature"));
     const std::vector<double> &T = field_vectors.at(field::temperature);
     for (unsigned int i = 0; i < property_vector.size(); ++i)
       property_vector[i] = A + B * T[i];
@@ -217,6 +223,9 @@ class ThermalConductivityPhaseChange : public ThermalConductivityModel
   double
   value(const std::map<field, double> &fields_value) override
   {
+    AssertThrow(fields_value.find(field::temperature) != fields_value.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "ThermalConductivityPhaseChange", "temperature"));
     double thermal_conductivity;
     // Thermal conductivity of solid phase
     if (fields_value.at(field::temperature) < T_solidus)
@@ -247,6 +256,9 @@ class ThermalConductivityPhaseChange : public ThermalConductivityModel
   vector_value(const std::map<field, std::vector<double>> &field_vectors,
                std::vector<double> &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "ThermalConductivityPhaseChange", "temperature"));
     const std::vector<double> &T = field_vectors.at(field::temperature);
     for (unsigned int i = 0; i < property_vector.size(); ++i)
       {
@@ -279,7 +291,12 @@ class ThermalConductivityPhaseChange : public ThermalConductivityModel
   jacobian(const std::map<field, double> &field_values, field id) override
   {
     if (id == field::temperature)
-      return this->numerical_jacobian(field_values, field::temperature);
+      {
+        AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                    PhysicialPropertyModelFieldUndefined(
+                      "ThermalConductivityPhaseChange", "temperature"));
+        return this->numerical_jacobian(field_values, field::temperature);
+      }
     else
       return 0;
   };
@@ -296,6 +313,9 @@ class ThermalConductivityPhaseChange : public ThermalConductivityModel
                   const field                                 id,
                   std::vector<double> &jacobian_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "ThermalConductivityPhaseChange", "temperature"));
     vector_numerical_jacobian(field_vectors, id, jacobian_vector);
   };
 
diff --git a/include/core/thermal_expansion_model.h b/include/core/thermal_expansion_model.h
index 87fccb2f5c..92eae2fd3d 100644
--- a/include/core/thermal_expansion_model.h
+++ b/include/core/thermal_expansion_model.h
@@ -134,6 +134,9 @@ class ThermalExpansionPhaseChange : public ThermalExpansionModel
   double
   value(const std::map<field, double> &fields_value) override
   {
+    AssertThrow(fields_value.find(field::temperature) != fields_value.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "ThermalExpansionPhaseChange", "temperature"));
     if (fields_value.at(field::temperature) > p_phase_change_params.T_liquidus)
       return p_phase_change_params.thermal_expansion_l;
     else
@@ -149,6 +152,9 @@ class ThermalExpansionPhaseChange : public ThermalExpansionModel
   vector_value(const std::map<field, std::vector<double>> &field_vectors,
                std::vector<double> &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "ThermalExpansionPhaseChange", "temperature"));
     const std::vector<double> &T = field_vectors.at(field::temperature);
     for (unsigned int i = 0; i < property_vector.size(); ++i)
       {
@@ -173,7 +179,12 @@ class ThermalExpansionPhaseChange : public ThermalExpansionModel
   jacobian(const std::map<field, double> &field_values, field id) override
   {
     if (id == field::temperature)
-      return this->numerical_jacobian(field_values, field::temperature);
+      {
+        AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                    PhysicialPropertyModelFieldUndefined(
+                      "ThermalExpansionPhaseChange", "temperature"));
+        return this->numerical_jacobian(field_values, field::temperature);
+      }
     else
       return 0;
   };
@@ -190,6 +201,9 @@ class ThermalExpansionPhaseChange : public ThermalExpansionModel
                   const field                                 id,
                   std::vector<double> &jacobian_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "ThermalExpansionPhaseChange", "temperature"));
     vector_numerical_jacobian(field_vectors, id, jacobian_vector);
   };
 
diff --git a/source/core/rheological_model.cc b/source/core/rheological_model.cc
index 0c8d50dd49..b78aa75ba0 100644
--- a/source/core/rheological_model.cc
+++ b/source/core/rheological_model.cc
@@ -39,6 +39,8 @@ Newtonian::value(const std::map<field, double> & /*field_values*/)
 double
 PowerLaw::value(const std::map<field, double> &field_values)
 {
+  AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
+              PhysicialPropertyModelFieldUndefined("PowerLaw", "shear_rate"));
   const double shear_rate_magnitude = field_values.at(field::shear_rate);
 
   return calculate_kinematic_viscosity(shear_rate_magnitude);
@@ -49,6 +51,8 @@ PowerLaw::vector_value(
   const std::map<field, std::vector<double>> &field_vectors,
   std::vector<double>                        &property_vector)
 {
+  AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+              PhysicialPropertyModelFieldUndefined("PowerLaw", "shear_rate"));
   const auto shear_rate_magnitude = field_vectors.at(field::shear_rate);
 
   for (unsigned int i = 0; i < shear_rate_magnitude.size(); ++i)
@@ -60,7 +64,12 @@ PowerLaw::jacobian(const std::map<field, double> &field_values, const field id)
 {
   const double shear_rate_magnitude = field_values.at(field::shear_rate);
   if (id == field::shear_rate)
-    return calculate_derivative(shear_rate_magnitude);
+    {
+      AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
+                  PhysicialPropertyModelFieldUndefined("PowerLaw",
+                                                       "shear_rate"));
+      return calculate_derivative(shear_rate_magnitude);
+    }
   else
     return 0;
 }
@@ -71,6 +80,8 @@ PowerLaw::vector_jacobian(
   const field                                 id,
   std::vector<double>                        &jacobian_vector)
 {
+  AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+              PhysicialPropertyModelFieldUndefined("PowerLaw", "shear_rate"));
   const auto shear_rate_magnitude = field_vectors.at(field::shear_rate);
 
   if (id == field::shear_rate)
@@ -83,6 +94,8 @@ PowerLaw::vector_jacobian(
 double
 Carreau::value(const std::map<field, double> &field_values)
 {
+  AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
+              PhysicialPropertyModelFieldUndefined("Carreau", "shear_rate"));
   const double shear_rate_magnitude = field_values.at(field::shear_rate);
 
   return calculate_kinematic_viscosity(shear_rate_magnitude);
@@ -92,6 +105,8 @@ void
 Carreau::vector_value(const std::map<field, std::vector<double>> &field_vectors,
                       std::vector<double> &property_vector)
 {
+  AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+              PhysicialPropertyModelFieldUndefined("Carreau", "shear_rate"));
   const auto shear_rate_magnitude = field_vectors.at(field::shear_rate);
 
   for (unsigned int i = 0; i < shear_rate_magnitude.size(); ++i)
@@ -105,7 +120,12 @@ double
 Carreau::jacobian(const std::map<field, double> &field_values, const field id)
 {
   if (id == field::shear_rate)
-    return this->numerical_jacobian(field_values, field::shear_rate);
+    {
+      AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
+                  PhysicialPropertyModelFieldUndefined("Carreau",
+                                                       "shear_rate"));
+      return this->numerical_jacobian(field_values, field::shear_rate);
+    }
   else
     return 0;
 }
@@ -117,9 +137,14 @@ Carreau::vector_jacobian(
   std::vector<double>                        &jacobian_vector)
 {
   if (id == field::shear_rate)
-    this->vector_numerical_jacobian(field_vectors,
-                                    field::shear_rate,
-                                    jacobian_vector);
+    {
+      AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+                  PhysicialPropertyModelFieldUndefined("Carreau",
+                                                       "shear_rate"));
+      this->vector_numerical_jacobian(field_vectors,
+                                      field::shear_rate,
+                                      jacobian_vector);
+    }
   else
     std::fill(jacobian_vector.begin(), jacobian_vector.end(), 0);
 }
@@ -129,6 +154,9 @@ Carreau::vector_jacobian(
 double
 PhaseChangeRheology::value(const std::map<field, double> &field_values)
 {
+  AssertThrow(field_values.find(field::temperature) != field_values.end(),
+              PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                   "temperature"));
   const double temperature = field_values.at(field::temperature);
 
   return kinematic_viscosity(temperature);
@@ -139,6 +167,9 @@ PhaseChangeRheology::vector_value(
   const std::map<field, std::vector<double>> &field_vectors,
   std::vector<double>                        &property_vector)
 {
+  AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+              PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                   "temperature"));
   const std::vector<double> &temperature_vec =
     field_vectors.at(field::temperature);
 
@@ -153,7 +184,12 @@ PhaseChangeRheology::jacobian(const std::map<field, double> &field_values,
                               const field                    id)
 {
   if (id == field::temperature)
-    return this->numerical_jacobian(field_values, field::temperature);
+    {
+      AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                  PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                       "temperature"));
+      return this->numerical_jacobian(field_values, field::temperature);
+    }
   else
     return 0;
 }
@@ -164,5 +200,8 @@ PhaseChangeRheology::vector_jacobian(
   const field                                 id,
   std::vector<double>                        &jacobian_vector)
 {
+  AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+              PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                   "temperature"));
   vector_numerical_jacobian(field_vectors, id, jacobian_vector);
 }