diff --git a/ProcessLib/TH2M/LocalAssemblerInterface.h b/ProcessLib/TH2M/LocalAssemblerInterface.h
index 4532720199b727f70ad67cd3351b00bc619f4de1..3c1cf69777ea88ef08dab9c834338190a29e0411 100644
--- a/ProcessLib/TH2M/LocalAssemblerInterface.h
+++ b/ProcessLib/TH2M/LocalAssemblerInterface.h
@@ -46,12 +46,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
         std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
         std::vector<double>& cache) const = 0;
 
-    virtual std::vector<double> const& getIntPtLiquidPressure(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
     virtual std::vector<double> const& getIntPtLiquidDensity(
         const double t,
         std::vector<GlobalVector*> const& x,
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 5449a3e40e1845ea95fb2e065690890c1951e77b..92101b0b9f8ab04e513e64f7da7d34edd1463514 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -37,9 +37,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
     : _process_data(process_data),
       _integration_method(integration_order),
       _element(e),
-      _is_axially_symmetric(is_axially_symmetric),
-      _liquid_pressure(
-          std::vector<double>(_integration_method.getNumberOfPoints()))
+      _is_axially_symmetric(is_axially_symmetric)
 {
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
@@ -318,15 +316,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
 
     assert(local_x.size() == matrix_size);
 
-    auto gas_pressure =
-        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
-            gas_pressure_size> const>(local_x.data() + gas_pressure_index,
-                                      gas_pressure_size);
-
-    auto capillary_pressure =
-        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
-            capillary_pressure_size> const>(
-            local_x.data() + capillary_pressure_index, capillary_pressure_size);
     updateConstitutiveVariables(
         Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()), t,
         std::numeric_limits<double>::quiet_NaN());
@@ -337,12 +326,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
     {
         auto& ip_data = _ip_data[ip];
         ip_data.pushBackState();
-
-        double const pGR = ip_data.N_p.dot(gas_pressure);
-        double const pCap = ip_data.N_p.dot(capillary_pressure);
-        double const pLR = pGR - pCap;
-
-        _liquid_pressure[ip] = pLR;
     }
 }
 
@@ -526,9 +509,7 @@ void TH2MLocalAssembler<
         auto& eps = ip.eps;
         auto const& sigma_eff = ip.sigma_eff;
 
-        double const pGR = Np.dot(gas_pressure);
         double const pCap = Np.dot(capillary_pressure);
-        double const pLR = pGR - pCap;
 
         GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
         GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
@@ -571,8 +552,6 @@ void TH2MLocalAssembler<
         // effective density
         auto const rho = phi_G * ip.rhoGR + phi_L * ip.rhoLR + phi_S * rho_SR;
 
-        _liquid_pressure[int_point] = pLR;
-
         // abbreviations
         const double rho_C_FR = s_G * ip.rhoCGR + s_L * ip.rhoCLR;
         const double rho_W_FR = s_G * ip.rhoWGR + s_L * ip.rhoWLR;
@@ -1099,6 +1078,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
     auto const capillary_pressure =
         local_x.template segment<capillary_pressure_size>(
             capillary_pressure_index);
+    auto const liquid_pressure = (gas_pressure - capillary_pressure).eval();
 
     NumLib::interpolateToHigherOrderNodes<
         ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
@@ -1110,6 +1090,11 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         DisplacementDim>(_element, _is_axially_symmetric, capillary_pressure,
                          *_process_data.capillary_pressure_interpolated);
 
+    NumLib::interpolateToHigherOrderNodes<
+        ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
+        DisplacementDim>(_element, _is_axially_symmetric, liquid_pressure,
+                         *_process_data.liquid_pressure_interpolated);
+
     auto const temperature =
         local_x.template segment<temperature_size>(temperature_index);
 
@@ -1127,17 +1112,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
 
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
-        auto& ip_data = _ip_data[ip];
-
-        auto const& Np = ip_data.N_p;
-
-        double const pGR = Np.dot(gas_pressure);
-        double const pCap = Np.dot(capillary_pressure);
-        double const pLR = pGR - pCap;
-
-        _liquid_pressure[ip] = pLR;
-
-        saturation_avg += ip_data.s_L;
+        saturation_avg += _ip_data[ip].s_L;
     }
     saturation_avg /= n_integration_points;
     (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index ece8445dc4a90e6493f262a64b5f2cf00f7682e4..aaa050e7f1c23804270c50afd03a635f4b8905e6 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -212,16 +212,6 @@ private:
     // be stored in some container to avoid defining one method for each
     // variable.
 
-    virtual std::vector<double> const& getIntPtLiquidPressure(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(!_liquid_pressure.empty());
-        return _liquid_pressure;
-    }
-
     virtual std::vector<double> const& getIntPtLiquidDensity(
         const double /*t*/,
         std::vector<GlobalVector*> const& /*x*/,
@@ -412,8 +402,6 @@ private:
         typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
         _secondary_data;
 
-    std::vector<double> _liquid_pressure;
-
     // The shape function of pressure has the same form with the shape function
     // of temperature
     static const int gas_pressure_index = 0;
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index b6165a7b6742b4b02071e3c8be957fbdce42381e..c7a5f37d1fc67d687fdbb691a180dbedc7d41b51 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -185,8 +185,6 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
                            &LocalAssemblerInterface::getIntPtSolidDensity);
     add_secondary_variable("liquid_density", 1,
                            &LocalAssemblerInterface::getIntPtLiquidDensity);
-    add_secondary_variable("liquid_pressure", 1,
-                           &LocalAssemblerInterface::getIntPtLiquidPressure);
     add_secondary_variable("mole_fraction_gas", 1,
                            &LocalAssemblerInterface::getIntPtMoleFractionGas);
     add_secondary_variable("mass_fraction_gas", 1,
@@ -238,6 +236,11 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
             const_cast<MeshLib::Mesh&>(mesh), "capillary_pressure_interpolated",
             MeshLib::MeshItemType::Node, 1);
 
+    _process_data.liquid_pressure_interpolated =
+        MeshLib::getOrCreateMeshProperty<double>(
+            const_cast<MeshLib::Mesh&>(mesh), "liquid_pressure_interpolated",
+            MeshLib::MeshItemType::Node, 1);
+
     _process_data.temperature_interpolated =
         MeshLib::getOrCreateMeshProperty<double>(
             const_cast<MeshLib::Mesh&>(mesh), "temperature_interpolated",
diff --git a/ProcessLib/TH2M/TH2MProcessData.h b/ProcessLib/TH2M/TH2MProcessData.h
index 59e991b556c3ae96d4696cedaff4a7d1803ed9ae..21222506a84ee68c7e0cbbe17fdef9ceaa2cb764 100644
--- a/ProcessLib/TH2M/TH2MProcessData.h
+++ b/ProcessLib/TH2M/TH2MProcessData.h
@@ -67,6 +67,7 @@ struct TH2MProcessData
     MeshLib::PropertyVector<double>* gas_pressure_interpolated = nullptr;
     MeshLib::PropertyVector<double>* capillary_pressure_interpolated = nullptr;
     MeshLib::PropertyVector<double>* temperature_interpolated = nullptr;
+    MeshLib::PropertyVector<double>* liquid_pressure_interpolated = nullptr;
 
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 };
diff --git a/Tests/Data/TH2M/M/M_2d_neumann/M_2d_neumann.prj b/Tests/Data/TH2M/M/M_2d_neumann/M_2d_neumann.prj
index 9a7f5f7be739713f8ca49f208db860a9281b7016..a0e822cc195b5947e04f6d5f917b5dfc25ac0991 100644
--- a/Tests/Data/TH2M/M/M_2d_neumann/M_2d_neumann.prj
+++ b/Tests/Data/TH2M/M/M_2d_neumann/M_2d_neumann.prj
@@ -48,7 +48,6 @@
                 <secondary_variable internal_name="velocity_liquid" output_name="velocity_liquid"/>
                 <secondary_variable internal_name="sigma" output_name="sigma"/>
                 <secondary_variable internal_name="epsilon" output_name="epsilon"/>
-                <secondary_variable internal_name="liquid_pressure" output_name="liquid_pressure"/>
                 <secondary_variable internal_name="liquid_density" output_name="liquid_density"/>
                 <secondary_variable internal_name="gas_density" output_name="gas_density"/>
                 <secondary_variable internal_name="porosity" output_name="porosity"/>
@@ -259,6 +258,7 @@
                 <variable>gas_pressure_interpolated</variable>
                 <variable>capillary_pressure</variable>
                 <variable>capillary_pressure_interpolated</variable>
+                <variable>liquid_pressure_interpolated</variable>
                 <variable>temperature</variable>
                 <variable>temperature_interpolated</variable>
                 <variable>displacement</variable>
@@ -266,7 +266,6 @@
                 <variable>epsilon</variable>
                 <variable>velocity_gas</variable>
                 <variable>velocity_liquid</variable>
-                <variable>liquid_pressure</variable>
                 <variable>liquid_density</variable>
                 <variable>gas_density</variable>
                 <variable>porosity</variable>
diff --git a/Tests/Data/TH2M/T/T_1d_dirichlet/T_1d_dirichlet.prj b/Tests/Data/TH2M/T/T_1d_dirichlet/T_1d_dirichlet.prj
index 8794bd52e038d69dbc88967d8edd5b128603ebb1..8d77e7522a64f11dadbb1581015b3435d51d9e43 100644
--- a/Tests/Data/TH2M/T/T_1d_dirichlet/T_1d_dirichlet.prj
+++ b/Tests/Data/TH2M/T/T_1d_dirichlet/T_1d_dirichlet.prj
@@ -48,7 +48,6 @@
                 <secondary_variable internal_name="velocity_liquid" output_name="velocity_liquid"/>
                 <secondary_variable internal_name="sigma" output_name="sigma"/>
                 <secondary_variable internal_name="epsilon" output_name="epsilon"/>
-                <secondary_variable internal_name="liquid_pressure" output_name="liquid_pressure"/>
                 <secondary_variable internal_name="liquid_density" output_name="liquid_density"/>
                 <secondary_variable internal_name="gas_density" output_name="gas_density"/>
                 <secondary_variable internal_name="porosity" output_name="porosity"/>
@@ -260,6 +259,7 @@
                 <variable>gas_pressure_interpolated</variable>
                 <variable>capillary_pressure</variable>
                 <variable>capillary_pressure_interpolated</variable>
+                <variable>liquid_pressure_interpolated</variable>
                 <variable>temperature</variable>
                 <variable>temperature_interpolated</variable>
                 <variable>displacement</variable>
@@ -267,7 +267,6 @@
                 <variable>epsilon</variable>
                 <variable>velocity_gas</variable>
                 <variable>velocity_liquid</variable>
-                <variable>liquid_pressure</variable>
                 <variable>liquid_density</variable>
                 <variable>gas_density</variable>
                 <variable>porosity</variable>
diff --git a/Tests/Data/TH2M/TH2/heatpipe_radial/heatpipe_radial.prj b/Tests/Data/TH2M/TH2/heatpipe_radial/heatpipe_radial.prj
index d854536ee1749bc49e879796661617e56a3d9a39..ad355f795cbeeaa88f4b18871c93bf191e9cf4d3 100644
--- a/Tests/Data/TH2M/TH2/heatpipe_radial/heatpipe_radial.prj
+++ b/Tests/Data/TH2M/TH2/heatpipe_radial/heatpipe_radial.prj
@@ -47,7 +47,6 @@
                 <secondary_variable internal_name="velocity_liquid" output_name="velocity_liquid"/>
                 <secondary_variable internal_name="sigma" output_name="sigma"/>
                 <secondary_variable internal_name="epsilon" output_name="epsilon"/>
-                <secondary_variable internal_name="liquid_pressure" output_name="liquid_pressure"/>
                 <secondary_variable internal_name="liquid_density" output_name="liquid_density"/>
                 <secondary_variable internal_name="gas_density" output_name="gas_density"/>
                 <secondary_variable internal_name="solid_density" output_name="solid_density"/>
@@ -334,6 +333,7 @@
                 <variable>gas_pressure_interpolated</variable>
                 <variable>capillary_pressure</variable>
                 <variable>capillary_pressure_interpolated</variable>
+                <variable>liquid_pressure_interpolated</variable>
                 <variable>temperature</variable>
                 <variable>temperature_interpolated</variable>
                 <variable>displacement</variable>
@@ -341,7 +341,6 @@
                 <variable>epsilon</variable>
                 <variable>velocity_gas</variable>
                 <variable>velocity_liquid</variable>
-                <variable>liquid_pressure</variable>
                 <variable>liquid_density</variable>
                 <variable>gas_density</variable>
                 <variable>solid_density</variable>
diff --git a/Tests/Data/TH2M/TH2/heatpipe_slab/heatpipe_slab.prj b/Tests/Data/TH2M/TH2/heatpipe_slab/heatpipe_slab.prj
index 0c43d5d1186f9b1d7b76cc00df036aec848d2525..f7f42d864afdec86661d3f9c9e0b23cafa829a5e 100644
--- a/Tests/Data/TH2M/TH2/heatpipe_slab/heatpipe_slab.prj
+++ b/Tests/Data/TH2M/TH2/heatpipe_slab/heatpipe_slab.prj
@@ -47,7 +47,6 @@
                 <secondary_variable internal_name="velocity_liquid" output_name="velocity_liquid"/>
                 <secondary_variable internal_name="sigma" output_name="sigma"/>
                 <secondary_variable internal_name="epsilon" output_name="epsilon"/>
-                <secondary_variable internal_name="liquid_pressure" output_name="liquid_pressure"/>
                 <secondary_variable internal_name="liquid_density" output_name="liquid_density"/>
                 <secondary_variable internal_name="gas_density" output_name="gas_density"/>
                 <secondary_variable internal_name="solid_density" output_name="solid_density"/>
@@ -344,6 +343,7 @@
                 <variable>gas_pressure_interpolated</variable>
                 <variable>capillary_pressure</variable>
                 <variable>capillary_pressure_interpolated</variable>
+                <variable>liquid_pressure_interpolated</variable>
                 <variable>temperature</variable>
                 <variable>temperature_interpolated</variable>
                 <variable>displacement</variable>
@@ -351,7 +351,6 @@
                 <variable>epsilon</variable>
                 <variable>velocity_gas</variable>
                 <variable>velocity_liquid</variable>
-                <variable>liquid_pressure</variable>
                 <variable>liquid_density</variable>
                 <variable>gas_density</variable>
                 <variable>solid_density</variable>
diff --git a/Tests/Data/TH2M/THM/slab/THM_1d_dirichlet.prj b/Tests/Data/TH2M/THM/slab/THM_1d_dirichlet.prj
index 76ce845f384e7554efde79f1a0feee690ac34cb8..ca15aa407a7167d14a87e12c66facccce6a46a18 100644
--- a/Tests/Data/TH2M/THM/slab/THM_1d_dirichlet.prj
+++ b/Tests/Data/TH2M/THM/slab/THM_1d_dirichlet.prj
@@ -46,7 +46,6 @@
                 <secondary_variable internal_name="velocity_liquid" output_name="velocity_liquid"/>
                 <secondary_variable internal_name="sigma" output_name="sigma"/>
                 <secondary_variable internal_name="epsilon" output_name="epsilon"/>
-                <secondary_variable internal_name="liquid_pressure" output_name="liquid_pressure"/>
                 <secondary_variable internal_name="liquid_density" output_name="liquid_density"/>
                 <secondary_variable internal_name="gas_density" output_name="gas_density"/>
                 <secondary_variable internal_name="porosity" output_name="porosity"/>
@@ -234,6 +233,7 @@
                 <variable>gas_pressure_interpolated</variable>
                 <variable>capillary_pressure</variable>
                 <variable>capillary_pressure_interpolated</variable>
+                <variable>liquid_pressure_interpolated</variable>
                 <variable>temperature</variable>
                 <variable>temperature_interpolated</variable>
                 <variable>displacement</variable>
@@ -241,7 +241,6 @@
                 <variable>epsilon</variable>
                 <variable>velocity_gas</variable>
                 <variable>velocity_liquid</variable>
-                <variable>liquid_pressure</variable>
                 <variable>liquid_density</variable>
                 <variable>gas_density</variable>
                 <variable>porosity</variable>
diff --git a/Tests/Data/TH2M/THM/sphere/point_heatsource.prj b/Tests/Data/TH2M/THM/sphere/point_heatsource.prj
index d6a05ed1d86130b358fa3fea0f28b7b5a555912a..e689725cca667c2d6fdebdc756ea82bcd727775c 100644
--- a/Tests/Data/TH2M/THM/sphere/point_heatsource.prj
+++ b/Tests/Data/TH2M/THM/sphere/point_heatsource.prj
@@ -47,7 +47,6 @@
                 <secondary_variable internal_name="velocity_liquid" output_name="velocity_liquid"/>
                 <secondary_variable internal_name="sigma" output_name="sigma"/>
                 <secondary_variable internal_name="epsilon" output_name="epsilon"/>
-                <secondary_variable internal_name="liquid_pressure" output_name="liquid_pressure"/>
                 <secondary_variable internal_name="liquid_density" output_name="liquid_density"/>
                 <secondary_variable internal_name="gas_density" output_name="gas_density"/>
                 <secondary_variable internal_name="porosity" output_name="porosity"/>
@@ -235,6 +234,7 @@
                 <variable>gas_pressure_interpolated</variable>
                 <variable>capillary_pressure</variable>
                 <variable>capillary_pressure_interpolated</variable>
+                <variable>liquid_pressure_interpolated</variable>
                 <variable>temperature</variable>
                 <variable>temperature_interpolated</variable>
                 <variable>displacement</variable>
@@ -242,7 +242,6 @@
                 <variable>epsilon</variable>
                 <variable>velocity_gas</variable>
                 <variable>velocity_liquid</variable>
-                <variable>liquid_pressure</variable>
                 <variable>liquid_density</variable>
                 <variable>gas_density</variable>
                 <variable>porosity</variable>
diff --git a/Tests/Data/TH2M/TH_unsaturated/heatpipe_radial_static_gas/heatpipe_radial_static_gas.prj b/Tests/Data/TH2M/TH_unsaturated/heatpipe_radial_static_gas/heatpipe_radial_static_gas.prj
index 77372cb987f2db91e33935512db6f8d92125d741..02ebf0ec183aed9d52323b7a2169157d1f7e63c0 100644
--- a/Tests/Data/TH2M/TH_unsaturated/heatpipe_radial_static_gas/heatpipe_radial_static_gas.prj
+++ b/Tests/Data/TH2M/TH_unsaturated/heatpipe_radial_static_gas/heatpipe_radial_static_gas.prj
@@ -47,7 +47,6 @@
                 <secondary_variable internal_name="velocity_liquid" output_name="velocity_liquid"/>
                 <secondary_variable internal_name="sigma" output_name="sigma"/>
                 <secondary_variable internal_name="epsilon" output_name="epsilon"/>
-                <secondary_variable internal_name="liquid_pressure" output_name="liquid_pressure"/>
                 <secondary_variable internal_name="liquid_density" output_name="liquid_density"/>
                 <secondary_variable internal_name="gas_density" output_name="gas_density"/>
                 <secondary_variable internal_name="solid_density" output_name="solid_density"/>
@@ -334,6 +333,7 @@
                 <variable>gas_pressure_interpolated</variable>
                 <variable>capillary_pressure</variable>
                 <variable>capillary_pressure_interpolated</variable>
+                <variable>liquid_pressure_interpolated</variable>
                 <variable>temperature</variable>
                 <variable>temperature_interpolated</variable>
                 <variable>displacement</variable>
@@ -341,7 +341,6 @@
                 <variable>epsilon</variable>
                 <variable>velocity_gas</variable>
                 <variable>velocity_liquid</variable>
-                <variable>liquid_pressure</variable>
                 <variable>liquid_density</variable>
                 <variable>gas_density</variable>
                 <variable>solid_density</variable>
diff --git a/Tests/Data/TH2M/TH_unsaturated/heatpipe_slab_static_gas/heatpipe_slab_static_gas.prj b/Tests/Data/TH2M/TH_unsaturated/heatpipe_slab_static_gas/heatpipe_slab_static_gas.prj
index ee054d49af4f5d0d5dfcf078630ce7d1b6864579..306707456a97398e8211e6daf2e5b57192683569 100644
--- a/Tests/Data/TH2M/TH_unsaturated/heatpipe_slab_static_gas/heatpipe_slab_static_gas.prj
+++ b/Tests/Data/TH2M/TH_unsaturated/heatpipe_slab_static_gas/heatpipe_slab_static_gas.prj
@@ -47,7 +47,6 @@
                 <secondary_variable internal_name="velocity_liquid" output_name="velocity_liquid"/>
                 <secondary_variable internal_name="sigma" output_name="sigma"/>
                 <secondary_variable internal_name="epsilon" output_name="epsilon"/>
-                <secondary_variable internal_name="liquid_pressure" output_name="liquid_pressure"/>
                 <secondary_variable internal_name="liquid_density" output_name="liquid_density"/>
                 <secondary_variable internal_name="gas_density" output_name="gas_density"/>
                 <secondary_variable internal_name="solid_density" output_name="solid_density"/>
@@ -344,6 +343,7 @@
                 <variable>gas_pressure_interpolated</variable>
                 <variable>capillary_pressure</variable>
                 <variable>capillary_pressure_interpolated</variable>
+                <variable>liquid_pressure_interpolated</variable>
                 <variable>temperature</variable>
                 <variable>temperature_interpolated</variable>
                 <variable>displacement</variable>
@@ -351,7 +351,6 @@
                 <variable>epsilon</variable>
                 <variable>velocity_gas</variable>
                 <variable>velocity_liquid</variable>
-                <variable>liquid_pressure</variable>
                 <variable>liquid_density</variable>
                 <variable>gas_density</variable>
                 <variable>solid_density</variable>