diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_explicit_hm_coupling_in_unsaturated_zone.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_explicit_hm_coupling_in_unsaturated_zone.md
new file mode 100644
index 0000000000000000000000000000000000000000..f3ea0d422a3fc7052e1e3626a03814c554db739c
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_explicit_hm_coupling_in_unsaturated_zone.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::RichardsMechanics::RichardsMechanicsProcessData::explicit_hm_coupling_in_unsaturated_zone
diff --git a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
index 7cfe1ce39fb5ccc4dc4ba207dc96a886be22d082..11b21d687ab0f18d41abfdc80ab02f41277ba374 100644
--- a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
@@ -173,13 +173,18 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
         //! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__mass_lumping}
         config.getConfigParameter<bool>("mass_lumping", false);
 
+    auto const explicit_hm_coupling_in_unsaturated_zone =
+        //! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__explicit_hm_coupling_in_unsaturated_zone}
+        config.getConfigParameter<bool>("explicit_hm_coupling_in_unsaturated_zone", false);
+
     RichardsMechanicsProcessData<DisplacementDim> process_data{
         materialIDs(mesh),
         std::move(media_map),
         std::move(solid_constitutive_relations),
         initial_stress,
         specific_body_force,
-        mass_lumping};
+        mass_lumping,
+        explicit_hm_coupling_in_unsaturated_zone};
 
     SecondaryVariableCollection secondary_variables;
 
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 1a86a5ebaa80d75bb087e5b70115028b27057336..cad65d3b19dd2af28c9c279169e117ab542f7587 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -954,8 +954,16 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         //
         // pressure equation, displacement part.
         //
-        Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
-                         identity2.transpose() * B * w;
+        if (_process_data.explicit_hm_coupling_in_unsaturated_zone)
+        {
+            Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR * alpha *
+                             identity2.transpose() * B * w;
+        }
+        else
+        {
+            Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
+                             identity2.transpose() * B * w;
+        }
 
         //
         // pressure equation, pressure part.
@@ -998,11 +1006,14 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
              specific_storage_a_S * dS_L_dp_cap) /
             dt * N_p * w;
 
-        local_Jac
-            .template block<pressure_size, pressure_size>(pressure_index,
-                                                          pressure_index)
-            .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
-                          identity2.transpose() * B * u_dot * N_p * w;
+        if (!_process_data.explicit_hm_coupling_in_unsaturated_zone)
+        {
+            local_Jac
+                .template block<pressure_size, pressure_size>(pressure_index,
+                                                              pressure_index)
+                .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
+                              identity2.transpose() * B * u_dot * N_p * w;
+        }
 
         double const dk_rel_dS_l =
             medium->property(MPL::PropertyType::relative_permeability)
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
index 1762e2f8b2301448aeebbad2c435cb83234af868..e42b6e1c3b95d3b7820b5b99d5a30112ba57ec1f 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
@@ -57,6 +57,10 @@ struct RichardsMechanicsProcessData
 
     bool const apply_mass_lumping;
 
+    /// If set, improves convergence of the global Newton method in unsaturated
+    /// regime.
+    bool const explicit_hm_coupling_in_unsaturated_zone;
+
     MeshLib::PropertyVector<double>* element_saturation = nullptr;
     MeshLib::PropertyVector<double>* element_porosity = nullptr;
     MeshLib::PropertyVector<double>* element_stresses = nullptr;