diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 154248a716b36a6b3576d6aab655ab79d156e17d..42d127276f9a65aaa8cf626ac7393123d3c2988d 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -1910,23 +1910,10 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                          temperature,
                          *this->process_data_.temperature_interpolated);
 
-    unsigned const n_integration_points =
-        this->integration_method_.getNumberOfPoints();
-
-    double saturation_avg = 0;
-
     ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
         this->solid_material_, *this->process_data_.phase_transition_model_};
 
     updateConstitutiveVariables(local_x, local_x_prev, t, dt, models);
-
-    for (unsigned ip = 0; ip < n_integration_points; ip++)
-    {
-        saturation_avg += this->current_states_[ip].S_L_data.S_L;
-    }
-    saturation_avg /= n_integration_points;
-    (*this->process_data_.element_saturation)[this->element_.getID()] =
-        saturation_avg;
 }
 
 }  // namespace TH2M
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index b3289ca9d146d002cf4841f474efc27f61ee5d8d..dc583602d1deb4f5410354e82a1077936b6a3c49 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -175,10 +175,6 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
             _process_data.solid_materials, local_assemblers_,
             _integration_point_writer, integration_order);
 
-    _process_data.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
-        const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
-        MeshLib::MeshItemType::Cell, 1);
-
     _process_data.gas_pressure_interpolated =
         MeshLib::getOrCreateMeshProperty<double>(
             const_cast<MeshLib::Mesh&>(mesh), "gas_pressure_interpolated",
@@ -320,6 +316,9 @@ void TH2MProcess<DisplacementDim>::computeSecondaryVariableConcrete(
         &LocalAssemblerInterface<DisplacementDim>::computeSecondaryVariable,
         local_assemblers_, pv.getActiveElementIDs(), getDOFTables(x.size()), t,
         dt, x, x_prev, process_id);
+
+    cell_average_data_.computeSecondaryVariable(DisplacementDim,
+                                                local_assemblers_);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/TH2M/TH2MProcessData.h b/ProcessLib/TH2M/TH2MProcessData.h
index 8aa138d9758cbc32bc70827c3eba720b67e2d77e..e23f42d8cff647772698b2a9a1237d88c4fc4c9a 100644
--- a/ProcessLib/TH2M/TH2MProcessData.h
+++ b/ProcessLib/TH2M/TH2MProcessData.h
@@ -52,7 +52,6 @@ struct TH2MProcessData
 
     const bool use_TaylorHood_elements;
 
-    MeshLib::PropertyVector<double>* element_saturation = nullptr;
     MeshLib::PropertyVector<double>* gas_pressure_interpolated = nullptr;
     MeshLib::PropertyVector<double>* capillary_pressure_interpolated = nullptr;
     MeshLib::PropertyVector<double>* temperature_interpolated = nullptr;