diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 91711e29547dd4daa8cd402bedad5c0504416944..d1515e5d2441256872a341544ed1c77ca604bf2e 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -51,29 +51,30 @@ std::ostream& operator<<(std::ostream& os,
 }
 
 template <typename Reactant>
-void setReactantAmount(Reactant& reactant,
-                       GlobalIndexType const& chemical_system_id,
-                       MaterialPropertyLib::Phase const& solid_phase,
-                       ParameterLib::SpatialPosition const& pos, double const t)
+void setReactantMolality(Reactant& reactant,
+                         GlobalIndexType const& chemical_system_id,
+                         MaterialPropertyLib::Phase const& solid_phase,
+                         ParameterLib::SpatialPosition const& pos,
+                         double const t)
 {
     auto const& solid_constituent = solid_phase.component(reactant.name);
 
-    auto const amount =
-        solid_constituent.property(MaterialPropertyLib::PropertyType::amount)
+    auto const molality =
+        solid_constituent.property(MaterialPropertyLib::PropertyType::molality)
             .template initialValue<double>(pos, t);
 
-    (*reactant.amount)[chemical_system_id] = amount;
+    (*reactant.molality)[chemical_system_id] = molality;
 }
 
 template <typename Reactant>
-static double averageReactantAmount(
+static double averageReactantMolality(
     Reactant const& reactant,
     std::vector<GlobalIndexType> const& chemical_system_indices)
 {
     double const sum = std::accumulate(
         chemical_system_indices.begin(), chemical_system_indices.end(), 0.0,
         [&](double const s, GlobalIndexType const id) {
-            return s + (*reactant.amount)[id];
+            return s + (*reactant.molality)[id];
         });
     return sum / chemical_system_indices.size();
 }
@@ -151,14 +152,14 @@ void PhreeqcIO::initializeChemicalSystemConcrete(
     auto const& solid_phase = medium->phase("Solid");
     for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
     {
-        setReactantAmount(kinetic_reactant, chemical_system_id, solid_phase,
-                          pos, t);
+        setReactantMolality(kinetic_reactant, chemical_system_id, solid_phase,
+                            pos, t);
     }
 
     for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
     {
-        setReactantAmount(equilibrium_reactant, chemical_system_id, solid_phase,
-                          pos, t);
+        setReactantMolality(equilibrium_reactant, chemical_system_id,
+                            solid_phase, pos, t);
     }
 }
 
@@ -530,7 +531,7 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
                         equilibrium_reactants.end(), compare_by_name,
                         "Could not find equilibrium reactant '" + item_name +
                             "'.");
-                    (*equilibrium_reactant.amount)[chemical_system_id] =
+                    (*equilibrium_reactant.molality)[chemical_system_id] =
                         accepted_items[item_id];
                     break;
                 }
@@ -541,7 +542,7 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
                         kinetic_reactants.begin(), kinetic_reactants.end(),
                         compare_by_name,
                         "Could not find kinetic reactant '" + item_name + "'.");
-                    (*kinetic_reactant.amount)[chemical_system_id] =
+                    (*kinetic_reactant.molality)[chemical_system_id] =
                         accepted_items[item_id];
                     break;
                 }
@@ -585,14 +586,15 @@ void PhreeqcIO::computeSecondaryVariable(
 {
     for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
     {
-        (*kinetic_reactant.amount_avg)[ele_id] =
-            averageReactantAmount(kinetic_reactant, chemical_system_indices);
+        (*kinetic_reactant.mesh_prop_molality)[ele_id] =
+            averageReactantMolality(kinetic_reactant, chemical_system_indices);
     }
 
     for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
     {
-        (*equilibrium_reactant.amount_avg)[ele_id] = averageReactantAmount(
-            equilibrium_reactant, chemical_system_indices);
+        (*equilibrium_reactant.mesh_prop_molality)[ele_id] =
+            averageReactantMolality(equilibrium_reactant,
+                                    chemical_system_indices);
     }
 }
 }  // namespace PhreeqcIOData
diff --git a/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp b/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
index 669c57d3b88816b9e901e3fbde076554adc0dc88..6a7286fb1c9012b5e6fd995076e727dd8aa12b7a 100644
--- a/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
+++ b/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
@@ -38,12 +38,12 @@ void ChemicalSystem::initialize(std::size_t const num_chemical_systems)
 
     for (auto& kinetic_reactant : kinetic_reactants)
     {
-        kinetic_reactant.amount->resize(num_chemical_systems);
+        kinetic_reactant.molality->resize(num_chemical_systems);
     }
 
     for (auto& equilibrium_reactant : equilibrium_reactants)
     {
-        equilibrium_reactant.amount->resize(num_chemical_systems);
+        equilibrium_reactant.molality->resize(num_chemical_systems);
     }
 }
 }  // namespace PhreeqcIOData
diff --git a/ChemistryLib/PhreeqcIOData/CreateEquilibriumReactants.cpp b/ChemistryLib/PhreeqcIOData/CreateEquilibriumReactants.cpp
index 198f4f59896f542c03487e0667c8e4629ddd6f1d..8455a6f23a95b11aeae478540dba07185a1b77b0 100644
--- a/ChemistryLib/PhreeqcIOData/CreateEquilibriumReactants.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateEquilibriumReactants.cpp
@@ -42,17 +42,15 @@ std::vector<EquilibriumReactant> createEquilibriumReactants(
             equilibrium_reactant_config.getConfigParameter<double>(
                 "saturation_index");
 
-        auto amount = MeshLib::getOrCreateMeshProperty<double>(
+        auto molality = MeshLib::getOrCreateMeshProperty<double>(
             mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
-        auto mesh_prop_amount = MeshLib::getOrCreateMeshProperty<double>(
+        auto mesh_prop_molality = MeshLib::getOrCreateMeshProperty<double>(
             mesh, name + "_avg", MeshLib::MeshItemType::Cell, 1);
-        mesh_prop_amount->resize(mesh.getNumberOfElements());
+        mesh_prop_molality->resize(mesh.getNumberOfElements());
 
-        equilibrium_reactants.emplace_back(std::move(name),
-                                           amount,
-                                           mesh_prop_amount,
-                                           saturation_index);
+        equilibrium_reactants.emplace_back(
+            std::move(name), molality, mesh_prop_molality, saturation_index);
     }
 
     return equilibrium_reactants;
diff --git a/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp b/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp
index a063b340b30e8ab12705f92284d56faebd64d41e..8ab72618334ff9349f37e8023baa4f29e1c42e16 100644
--- a/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp
@@ -50,12 +50,12 @@ std::vector<KineticReactant> createKineticReactants(
             //! \ogs_file_param{prj__chemical_system__kinetic_reactants__kinetic_reactant__fix_amount}
             reactant_config.getConfigParameter<bool>("fix_amount", false);
 
-        auto amount = MeshLib::getOrCreateMeshProperty<double>(
+        auto molality = MeshLib::getOrCreateMeshProperty<double>(
             mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
-        auto mesh_prop_amount = MeshLib::getOrCreateMeshProperty<double>(
+        auto mesh_prop_molality = MeshLib::getOrCreateMeshProperty<double>(
             mesh, name + "_avg", MeshLib::MeshItemType::Cell, 1);
-        mesh_prop_amount->resize(mesh.getNumberOfElements());
+        mesh_prop_molality->resize(mesh.getNumberOfElements());
 
         if (chemical_formula.empty() && fix_amount)
         {
@@ -66,8 +66,8 @@ std::vector<KineticReactant> createKineticReactants(
 
         kinetic_reactants.emplace_back(std::move(name),
                                        std::move(chemical_formula),
-                                       amount,
-                                       mesh_prop_amount,
+                                       molality,
+                                       mesh_prop_molality,
                                        std::move(parameters),
                                        fix_amount);
     }
diff --git a/ChemistryLib/PhreeqcIOData/EquilibriumReactant.cpp b/ChemistryLib/PhreeqcIOData/EquilibriumReactant.cpp
index f6ff1020f85de0c89aa85888067b471f80b43146..89c88653a11fa822b43f729a0de92614386402c9 100644
--- a/ChemistryLib/PhreeqcIOData/EquilibriumReactant.cpp
+++ b/ChemistryLib/PhreeqcIOData/EquilibriumReactant.cpp
@@ -19,8 +19,8 @@ namespace PhreeqcIOData
 void EquilibriumReactant::print(std::ostream& os,
                                 std::size_t const global_id) const
 {
-    os << name << " " << saturation_index << " "
-       << (*amount)[global_id] << "\n";
+    os << name << " " << saturation_index << " " << (*molality)[global_id]
+       << "\n";
 }
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/EquilibriumReactant.h b/ChemistryLib/PhreeqcIOData/EquilibriumReactant.h
index 3ef18ab9e9c82bbff80790eb81f7c5b812d07637..c502f756d0577d5802de6f31b6cdf7a6c4b7af0e 100644
--- a/ChemistryLib/PhreeqcIOData/EquilibriumReactant.h
+++ b/ChemistryLib/PhreeqcIOData/EquilibriumReactant.h
@@ -29,12 +29,12 @@ namespace PhreeqcIOData
 struct EquilibriumReactant
 {
     EquilibriumReactant(std::string name_,
-                        MeshLib::PropertyVector<double>* amount_,
-                        MeshLib::PropertyVector<double>* amount_avg_,
+                        MeshLib::PropertyVector<double>* molality_,
+                        MeshLib::PropertyVector<double>* mesh_prop_molality_,
                         double saturation_index_)
         : name(std::move(name_)),
-          amount(amount_),
-          amount_avg(amount_avg_),
+          molality(molality_),
+          mesh_prop_molality(mesh_prop_molality_),
           saturation_index(saturation_index_)
     {
     }
@@ -42,8 +42,8 @@ struct EquilibriumReactant
     void print(std::ostream& os, std::size_t const global_id) const;
 
     std::string const name;
-    MeshLib::PropertyVector<double>* amount;
-    MeshLib::PropertyVector<double>* amount_avg;
+    MeshLib::PropertyVector<double>* molality;
+    MeshLib::PropertyVector<double>* mesh_prop_molality;
     double const saturation_index;
     static const ItemType item_type = ItemType::EquilibriumReactant;
 };
diff --git a/ChemistryLib/PhreeqcIOData/KineticReactant.cpp b/ChemistryLib/PhreeqcIOData/KineticReactant.cpp
index 18db0961fe7272fb3b01a927c9edb5923aea9452..ac8b6fe83d81b87c921c983f05f1107a08be4ebd 100644
--- a/ChemistryLib/PhreeqcIOData/KineticReactant.cpp
+++ b/ChemistryLib/PhreeqcIOData/KineticReactant.cpp
@@ -26,7 +26,7 @@ void KineticReactant::print(std::ostream& os,
         os << "-formula " << chemical_formula << "\n";
     }
 
-    os << "-m  " << (*amount)[global_id] << "\n";
+    os << "-m  " << (*molality)[global_id] << "\n";
 
     if (!parameters.empty())
     {
diff --git a/ChemistryLib/PhreeqcIOData/KineticReactant.h b/ChemistryLib/PhreeqcIOData/KineticReactant.h
index 40c9af5b6f4f7fe9c16e7c03bafeecfae20aeaf5..b104206f71f66703be5924cad1d737a02ce45b54 100644
--- a/ChemistryLib/PhreeqcIOData/KineticReactant.h
+++ b/ChemistryLib/PhreeqcIOData/KineticReactant.h
@@ -26,14 +26,14 @@ struct KineticReactant
 {
     KineticReactant(std::string name_,
                     std::string chemical_formula_,
-                    MeshLib::PropertyVector<double>* amount_,
-                    MeshLib::PropertyVector<double>* amount_avg_,
+                    MeshLib::PropertyVector<double>* molality_,
+                    MeshLib::PropertyVector<double>* mesh_prop_molality_,
                     std::vector<double>&& parameters_,
                     bool const fix_amount_)
         : name(std::move(name_)),
           chemical_formula(std::move(chemical_formula_)),
-          amount(amount_),
-          amount_avg(amount_avg_),
+          molality(molality_),
+          mesh_prop_molality(mesh_prop_molality_),
           parameters(std::move(parameters_)),
           fix_amount(fix_amount_)
     {
@@ -43,8 +43,8 @@ struct KineticReactant
 
     std::string const name;
     std::string const chemical_formula;
-    MeshLib::PropertyVector<double>* amount;
-    MeshLib::PropertyVector<double>* amount_avg;
+    MeshLib::PropertyVector<double>* molality;
+    MeshLib::PropertyVector<double>* mesh_prop_molality;
     std::vector<double> const parameters;
     bool const fix_amount;
     static const ItemType item_type = ItemType::KineticReactant;
diff --git a/MaterialLib/MPL/PropertyType.h b/MaterialLib/MPL/PropertyType.h
index 55a7adb99c0793ef778d0367ab7337d8e05b5e96..554e7a52bb131f56752fe4667b407e6ef0228dc0 100644
--- a/MaterialLib/MPL/PropertyType.h
+++ b/MaterialLib/MPL/PropertyType.h
@@ -35,7 +35,6 @@ namespace MaterialPropertyLib
 enum PropertyType : int
 {
     acentric_factor,
-    amount,
     binary_interaction_coefficient,
     biot_coefficient,
     bishops_effective_stress,
@@ -58,6 +57,7 @@ enum PropertyType : int
     heat_capacity,
     /// used to compute the hydrodynamic dispersion tensor.
     longitudinal_dispersivity,
+    molality,
     molar_mass,
     mole_fraction,
     /// used to compute the hydrodynamic dispersion tensor.
@@ -182,6 +182,11 @@ inline PropertyType convertStringToProperty(std::string const& inString)
     {
         return PropertyType::longitudinal_dispersivity;
     }
+    // TODO (renchao): add property "volume fraction"
+    if (boost::iequals(inString, "molality"))
+    {
+        return PropertyType::molality;
+    }
     if (boost::iequals(inString, "molar_mass"))
     {
         return PropertyType::molar_mass;
@@ -290,11 +295,6 @@ inline PropertyType convertStringToProperty(std::string const& inString)
     {
         return PropertyType::viscosity;
     }
-    // TODO (renchao): replace property "amount" with volume fraction
-    if (boost::iequals(inString, "amount"))
-    {
-        return PropertyType::amount;
-    }
 
     OGS_FATAL(
         "The property name '{:s}' does not correspond to any known property",
@@ -306,7 +306,6 @@ inline PropertyType convertStringToProperty(std::string const& inString)
 
 static const std::array<std::string, PropertyType::number_of_properties>
     property_enum_to_string{{"acentric_factor",
-                             "amount",
                              "binary_interaction_coefficient",
                              "biot_coefficient",
                              "bishops_effective_stress",
@@ -327,6 +326,7 @@ static const std::array<std::string, PropertyType::number_of_properties>
                              "fredlund_parameters",
                              "heat_capacity",
                              "longitudinal_dispersivity",
+                             "molality",
                              "molar_mass",
                              "mole_fraction",
                              "molecular_diffusion",