diff --git a/ProcessLib/TH2M/ConstitutiveVariables.h b/ProcessLib/TH2M/ConstitutiveVariables.h
index 6a06e5a598c66a5c0085bc9964078b8e18055095..e7aaa38995b69015b8d95cac5e092d0700493ea4 100644
--- a/ProcessLib/TH2M/ConstitutiveVariables.h
+++ b/ProcessLib/TH2M/ConstitutiveVariables.h
@@ -47,6 +47,12 @@ struct ConstitutiveVariables
     DisplacementDimMatrix dfW_4_LWpC_d_dp_cap;
     DisplacementDimMatrix dfW_4_LWpC_d_dT;
     DisplacementDimMatrix dfC_4_LCpG_dT;
+    DisplacementDimMatrix dfC_4_LCpC_a_dp_GR;
+    DisplacementDimMatrix dfC_4_LCpC_a_dp_cap;
+    DisplacementDimMatrix dfC_4_LCpC_a_dT;
+    DisplacementDimMatrix dfC_4_LCpC_d_dp_GR;
+    DisplacementDimMatrix dfC_4_LCpC_d_dp_cap;
+    DisplacementDimMatrix dfC_4_LCpC_d_dT;
     DisplacementDimMatrix dadvection_C_dp_GR;
     DisplacementDimMatrix dadvection_C_dp_cap;
     DisplacementDimMatrix dk_over_mu_G_dp_cap;
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index fe85628cce329b62571866d935b1d7fc616f5912..b0f05849d61f772dc332f5b20d1b2aa6d68ffa9c 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -698,6 +698,23 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
         ip_cv.dfW_4_LWpC_d_dT =
             Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
+
+        ip_cv.dfC_4_LCpC_a_dp_GR = c.drho_C_LR_dp_GR * k_over_mu_L
+            //+ rhoCLR * (dk_over_mu_L_dp_GR = 0)
+            ;
+        ip_cv.dfC_4_LCpC_a_dp_cap = -c.drho_C_LR_dp_LR * k_over_mu_L +
+                                    ip_data.rhoCLR * ip_cv.dk_over_mu_L_dp_cap;
+        ip_cv.dfC_4_LCpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
+            //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
+            ;
+
+        // TODO (naumov) for dxmW*/d* != 0
+        ip_cv.dfC_4_LCpC_d_dp_GR =
+            Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
+        ip_cv.dfC_4_LCpC_d_dp_cap =
+            Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
+        ip_cv.dfC_4_LCpC_d_dT =
+            Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
     }
 
     return ip_constitutive_variables;
@@ -1492,6 +1509,32 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         LCpC.noalias() -=
             gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w;
 
+        /* TODO (naumov) This part is not tested by any of the current ctests.
+        // d (fC_4_LCpC * grad p_cap)/d p_GR
+        local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
+            gradNpT *
+            (ip_cv.dfC_4_LCpC_a_dp_GR
+             // + ip_cv.dfC_4_LCpC_d_dp_GR TODO (naumov)
+             ) *
+            gradpCap * Np * w;
+        // d (fC_4_LCpC * grad p_cap)/d p_cap
+        local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
+            gradNpT *
+            (ip_cv.dfC_4_LCpC_a_dp_cap
+             // + ip_cv.dfC_4_LCpC_d_dp_cap TODO (naumov)
+             ) *
+            gradpCap * Np * w;
+
+        local_Jac
+            .template block<C_size, temperature_size>(C_index,
+                                                      temperature_index)
+            .noalias() += gradNpT *
+                          (ip_cv.dfC_4_LCpC_a_dT
+                           // + ip_cv.dfC_4_LCpC_d_dT TODO (naumov)
+                           ) *
+                          gradpCap * Np * w;
+        */
+
         LCT.noalias() += gradNpT * diffusion_C_T * gradNp * w;
 
         // fC_1