diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp
index 1c7ba200bd..df99c5df9d 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp
@@ -993,7 +993,7 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh,
   } );
 
   // check if comp fractions need to be corrected to avoid zero diags etc
-  if( m_useZFormulation )
+  if( m_useZFormulation && m_allowCompDensChopping )
   {
     DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
     chopNegativeCompFractions( domain );
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
index e5636f6920..cc3062f171 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
@@ -672,8 +672,6 @@ real64 CompositionalMultiphaseFVM::scalingForSystemSolutionZFormulation( DomainP
                                                                          DofManager const & dofManager,
                                                                          arrayView1d< real64 const > const & localSolution )
 {
-  GEOS_MARK_FUNCTION;
-
   string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
   real64 scalingFactor = 1.0;
   real64 maxDeltaPres = 0.0, maxDeltaCompFrac = 0.0, maxDeltaTemp = 0.0;
@@ -954,8 +952,6 @@ void CompositionalMultiphaseFVM::applySystemSolution( DofManager const & dofMana
   // if component density chopping is allowed, some component densities may be negative after the update
   // these negative component densities are set to zero in this function
 
-  // TODO: implement chopNegativeZ
-
   if( m_allowCompDensChopping )
   {
     if( m_useZFormulation )
@@ -968,15 +964,8 @@ void CompositionalMultiphaseFVM::applySystemSolution( DofManager const & dofMana
                                                                MeshLevel & mesh,
                                                                arrayView1d< string const > const & regionNames )
   {
-    std::vector< string > fields{ fields::flow::pressure::key() };
-    if( m_useZFormulation )
-    {
-      fields.emplace_back( fields::flow::globalCompFraction::key() );
-    }
-    else
-    {
-      fields.emplace_back( fields::flow::globalCompDensity::key() );
-    }
+    std::vector< string > fields{ fields::flow::pressure::key(),
+                                  m_useZFormulation ? fields::flow::globalCompFraction::key() : fields::flow::globalCompDensity::key() };
     if( m_isThermal )
     {
       fields.emplace_back( fields::flow::temperature::key() );
@@ -1001,22 +990,14 @@ void CompositionalMultiphaseFVM::updatePhaseMobility( ObjectManagerBase & dataGr
 
   if( m_useZFormulation )
   {
-    if( m_isThermal )
-    {
-      // For now: isothermal only
-      GEOS_ERROR_IF( m_isThermal, GEOS_FMT(
-                       "{}: Z Formulation is currently not available for thermal simulations", getDataContext() ) );
-    }
-    else
-    {
-      isothermalCompositionalMultiphaseFVMKernels::
-        PhaseMobilityZFormulationKernelFactory::
-        createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
-                                                   m_numPhases,
-                                                   dataGroup,
-                                                   fluid,
-                                                   relperm );
-    }
+    // For now: isothermal only
+    isothermalCompositionalMultiphaseFVMKernels::
+      PhaseMobilityZFormulationKernelFactory::
+      createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                 m_numPhases,
+                                                 dataGroup,
+                                                 fluid,
+                                                 relperm );
   }
   else
   {
@@ -1261,8 +1242,6 @@ void CompositionalMultiphaseFVM::applyFaceDirichletBC( real64 const time_n,
 
       if( m_useZFormulation )
       {
-        GEOS_ERROR_IF( m_isThermal, GEOS_FMT(
-                         "{}: Z Formulation is currently not available for thermal simulations", getDataContext() ) );
         isothermalCompositionalMultiphaseFVMKernels::
           DirichletFluxComputeZFormulationKernelFactory::
           createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp
index 2798752e69..f1437bd45f 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp
@@ -109,7 +109,7 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase
                             DofManager const & dofManager,
                             arrayView1d< real64 const > const & localSolution ) override;
 
-  virtual real64
+  real64
   scalingForSystemSolutionZFormulation( DomainPartition & domain,
                                         DofManager const & dofManager,
                                         arrayView1d< real64 const > const & localSolution );
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp
index 909997e7d6..922b8fbf01 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp
@@ -509,17 +509,6 @@ real64 CompositionalMultiphaseHybridFVM::scalingForSystemSolution( DomainPartiti
   return LvArray::math::max( MpiWrapper::min( scalingFactor ), m_minScalingFactor );
 }
 
-real64 CompositionalMultiphaseHybridFVM::scalingForSystemSolutionZFormulation( DomainPartition & domain,
-                                                                               DofManager const & dofManager,
-                                                                               arrayView1d< real64 const > const & localSolution )
-{
-  GEOS_MARK_FUNCTION;
-  // z formulation not implemented
-  GEOS_UNUSED_VAR( domain, dofManager, localSolution );
-  GEOS_ERROR( "Z formulation not yet available for this flow solver" );
-  return 0;
-}
-
 bool CompositionalMultiphaseHybridFVM::checkSystemSolution( DomainPartition & domain,
                                                             DofManager const & dofManager,
                                                             arrayView1d< real64 const > const & localSolution,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp
index 5fa78d20ce..d133610a7b 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp
@@ -111,11 +111,6 @@ class CompositionalMultiphaseHybridFVM : public CompositionalMultiphaseBase
                             DofManager const & dofManager,
                             arrayView1d< real64 const > const & localSolution ) override;
 
-  virtual real64
-  scalingForSystemSolutionZFormulation( DomainPartition & domain,
-                                        DofManager const & dofManager,
-                                        arrayView1d< real64 const > const & localSolution );
-
   virtual bool
   checkSystemSolution( DomainPartition & domain,
                        DofManager const & dofManager,