diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp
index 66c7967703..df99c5df9d 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp
@@ -54,6 +54,8 @@
 #include "physicsSolvers/fluidFlow/kernels/compositional/HydrostaticPressureKernel.hpp"
 #include "physicsSolvers/fluidFlow/kernels/compositional/StatisticsKernel.hpp"
 #include "physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/zFormulation/AccumulationZFormulationKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseVolumeFractionZFormulationKernel.hpp"
 
 #if defined( __INTEL_COMPILER )
 #pragma GCC optimize "O0"
@@ -79,7 +81,8 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name,
   m_useTotalMassEquation( 1 ),
   m_useSimpleAccumulation( 1 ),
   m_useNewGravity( 0 ),
-  m_minCompDens( isothermalCompositionalMultiphaseBaseKernels::minDensForDivision )
+  m_minCompDens( isothermalCompositionalMultiphaseBaseKernels::minDensForDivision ),
+  m_minCompFrac( isothermalCompositionalMultiphaseBaseKernels::minCompFracForDivision )
 {
 //START_SPHINX_INCLUDE_00
   this->registerWrapper( viewKeyStruct::inputTemperatureString(), &m_inputTemperature ).
@@ -93,6 +96,10 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name,
     setInputFlag( InputFlags::OPTIONAL ).
     setDescription( GEOS_FMT( "Use mass formulation instead of molar. Warning : Affects {} rates units.",
                               SourceFluxBoundaryCondition::catalogName() ) );
+  this->registerWrapper( viewKeyStruct::useZFormulationFlagString(), &m_useZFormulation ).
+    setApplyDefaultValue( 0 ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setDescription( "Use overall composition (Z) formulation instead of component density." );
 
   this->registerWrapper( viewKeyStruct::solutionChangeScalingFactorString(), &m_solutionChangeScalingFactor ).
     setSizedFromParent( 0 ).
@@ -119,6 +126,11 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name,
     setInputFlag( InputFlags::OPTIONAL ).
     setApplyDefaultValue( LvArray::NumericLimits< real64 >::max ). // disabled by default
     setDescription( "Target (relative) change in component density in a time step" );
+  this->registerWrapper( viewKeyStruct::targetCompFracChangeString(), &m_targetCompFracChange ).
+    setSizedFromParent( 0 ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setApplyDefaultValue( LvArray::NumericLimits< real64 >::max ). // disabled by default
+    setDescription( "Target change in component fraction in a time step" );
 
   this->registerWrapper( viewKeyStruct::maxCompFracChangeString(), &m_maxCompFracChange ).
     setSizedFromParent( 0 ).
@@ -177,6 +189,12 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name,
     setApplyDefaultValue( isothermalCompositionalMultiphaseBaseKernels::minDensForDivision ).
     setDescription( "Minimum allowed global component density" );
 
+  this->registerWrapper( viewKeyStruct::minCompFracString(), &m_minCompFrac ).
+    setSizedFromParent( 0 ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setApplyDefaultValue( isothermalCompositionalMultiphaseBaseKernels::minCompFracForDivision ).
+    setDescription( "Minimum allowed global component fraction" );
+
   this->registerWrapper( viewKeyStruct::maxSequentialCompDensChangeString(), &m_maxSequentialCompDensChange ).
     setSizedFromParent( 0 ).
     setInputFlag( InputFlags::OPTIONAL ).
@@ -219,8 +237,11 @@ void CompositionalMultiphaseBase::postInputInitialization()
                         getWrapperDataContext( viewKeyStruct::targetPhaseVolFracChangeString() ) <<
                         ": The target change in phase volume fraction in a time step must be larger than 0.0" );
   GEOS_ERROR_IF_LE_MSG( m_targetRelativeCompDensChange, 0.0,
-                        getWrapperDataContext( viewKeyStruct::targetPhaseVolFracChangeString() ) <<
+                        getWrapperDataContext( viewKeyStruct::targetRelativeCompDensChangeString() ) <<
                         ": The target change in component density in a time step must be larger than 0.0" );
+  GEOS_ERROR_IF_LE_MSG( m_targetCompFracChange, 0.0,
+                        getWrapperDataContext( viewKeyStruct::targetCompFracChangeString() ) <<
+                        ": The target change in component fraction in a time step must be larger than 0.0" );
 
   GEOS_ERROR_IF_LT_MSG( m_solutionChangeScalingFactor, 0.0,
                         getWrapperDataContext( viewKeyStruct::solutionChangeScalingFactorString() ) <<
@@ -235,11 +256,31 @@ void CompositionalMultiphaseBase::postInputInitialization()
                         getWrapperDataContext( viewKeyStruct::minScalingFactorString() ) <<
                         ": The minumum scaling factor must be smaller or equal to 1.0" );
 
-  if( m_isThermal && m_useSimpleAccumulation == 1 ) // useSimpleAccumulation is not yet compatible with thermal
+  if( m_isThermal && m_useSimpleAccumulation ) // useSimpleAccumulation is not yet compatible with thermal
   {
-    GEOS_LOG_RANK_0( "'useSimpleAccumulation' is not yet implemented for thermal simulation. Switched to phase sum accumulation." );
+    GEOS_LOG_RANK_0( GEOS_FMT( "{}: '{}' is not yet implemented for thermal simulation. Switched to phase sum accumulation.",
+                               getDataContext(), viewKeyStruct::useSimpleAccumulationString() ) );
     m_useSimpleAccumulation = 0;
   }
+
+  if( m_useZFormulation )
+  {
+    if( m_isThermal ) // useZFormulation is not yet compatible with thermal
+    {
+      GEOS_ERROR( GEOS_FMT( "{}: '{}' is currently not available for thermal simulations",
+                            getDataContext(), viewKeyStruct::useZFormulationFlagString() ) );
+    }
+    if( m_hasDiffusion || m_hasDispersion )
+    {
+      GEOS_ERROR( GEOS_FMT( "{}: {} is currently not available for diffusion or dispersion",
+                            getDataContext(), viewKeyStruct::useZFormulationFlagString() ) );
+    }
+    if( m_isJumpStabilized ) // useZFormulation is not yet compatible with pressure stabilization
+    {
+      GEOS_ERROR( GEOS_FMT( "{}: pressure stabilization is not yet supported by {}",
+                            getDataContext(), viewKeyStruct::useZFormulationFlagString() ) );
+    }
+  }
 }
 
 void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies )
@@ -300,8 +341,21 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies )
     m_isThermal = referenceFluid.isThermal();
   }
 
-  // n_c components + one pressure ( + one temperature if needed )
-  m_numDofPerCell = m_isThermal ? m_numComponents + 2 : m_numComponents + 1;
+
+  if( m_useZFormulation )
+  {
+    // Z formulation - component densities and pressure are primary unknowns
+    // (n_c-1) overall compositions + one pressure
+    // Testing: let's have sum_c zc = 1 as explicit equation for now
+    m_numDofPerCell = m_numComponents + 1;
+  }
+  else
+  {
+    // default formulation - component densities and pressure are primary unknowns
+    // n_c components + one pressure ( + one temperature if needed )
+    m_numDofPerCell = m_isThermal ? m_numComponents + 2 : m_numComponents + 1;
+  }
+
 
   // 2. Register and resize all fields as necessary
   forDiscretizationOnMeshTargets( meshBodies, [&]( string const &,
@@ -374,28 +428,43 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies )
 
       subRegion.registerField< pressureScalingFactor >( getName() );
       subRegion.registerField< temperatureScalingFactor >( getName() );
-      subRegion.registerField< globalCompDensityScalingFactor >( getName() );
 
       // The resizing of the arrays needs to happen here, before the call to initializePreSubGroups,
       // to make sure that the dimensions are properly set before the timeHistoryOutput starts its initialization.
-
-      subRegion.registerField< globalCompDensity >( getName() ).
+      subRegion.registerField< globalCompFraction >( getName() ).
         setDimLabels( 1, fluid.componentNames() ).
         reference().resizeDimension< 1 >( m_numComponents );
-      subRegion.registerField< globalCompDensity_n >( getName() ).
-        reference().resizeDimension< 1 >( m_numComponents );
-      if( m_isFixedStressPoromechanicsUpdate )
+      if( m_useZFormulation )
       {
-        subRegion.registerField< globalCompDensity_k >( getName() ).
+        subRegion.registerField< globalCompFraction_n >( getName() ).
           setDimLabels( 1, fluid.componentNames() ).
           reference().resizeDimension< 1 >( m_numComponents );
+        // may be needed later for sequential poromechanics implementation
+        //if( m_isFixedStressPoromechanicsUpdate )
+        //{
+        //  subRegion.registerField< globalCompFraction_k >( getName() ).
+        //    setDimLabels( 1, fluid.componentNames() ).
+        //    reference().resizeDimension< 1 >( m_numComponents );
+        //}
+        subRegion.registerField< globalCompFractionScalingFactor >( getName() );
+      }
+      else
+      {
+        subRegion.registerField< globalCompDensity >( getName() ).
+          setDimLabels( 1, fluid.componentNames() ).
+          reference().resizeDimension< 1 >( m_numComponents );
+        subRegion.registerField< globalCompDensity_n >( getName() ).
+          reference().resizeDimension< 1 >( m_numComponents );
+        if( m_isFixedStressPoromechanicsUpdate )
+        {
+          subRegion.registerField< globalCompDensity_k >( getName() ).
+            setDimLabels( 1, fluid.componentNames() ).
+            reference().resizeDimension< 1 >( m_numComponents );
+        }
+        subRegion.registerField< globalCompDensityScalingFactor >( getName() );
+        subRegion.registerField< dGlobalCompFraction_dGlobalCompDensity >( getName() ).
+          reference().resizeDimension< 1, 2 >( m_numComponents, m_numComponents );
       }
-
-      subRegion.registerField< globalCompFraction >( getName() ).
-        setDimLabels( 1, fluid.componentNames() ).
-        reference().resizeDimension< 1 >( m_numComponents );
-      subRegion.registerField< dGlobalCompFraction_dGlobalCompDensity >( getName() ).
-        reference().resizeDimension< 1, 2 >( m_numComponents, m_numComponents );
 
       subRegion.registerField< phaseVolumeFraction >( getName() ).
         setDimLabels( 1, fluid.phaseNames() ).
@@ -658,22 +727,37 @@ real64 CompositionalMultiphaseBase::updatePhaseVolumeFraction( ObjectManagerBase
   string const & fluidName = dataGroup.getReference< string >( viewKeyStruct::fluidNamesString() );
   MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( dataGroup, fluidName );
 
-  real64 maxDeltaPhaseVolFrac  =
-    m_isThermal ?
-    thermalCompositionalMultiphaseBaseKernels::
-      PhaseVolumeFractionKernelFactory::
-      createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
-                                                 m_numPhases,
-                                                 dataGroup,
-                                                 fluid )
-:    isothermalCompositionalMultiphaseBaseKernels::
-      PhaseVolumeFractionKernelFactory::
-      createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
-                                                 m_numPhases,
-                                                 dataGroup,
-                                                 fluid );
-
-  return maxDeltaPhaseVolFrac;
+  if( m_useZFormulation )
+  {
+    // isothermal for now
+    return isothermalCompositionalMultiphaseBaseKernels::
+             PhaseVolumeFractionZFormulationKernelFactory::
+             createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                        m_numPhases,
+                                                        dataGroup,
+                                                        fluid );
+  }
+  else
+  {
+    if( m_isThermal )
+    {
+      return thermalCompositionalMultiphaseBaseKernels::
+               PhaseVolumeFractionKernelFactory::
+               createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                          m_numPhases,
+                                                          dataGroup,
+                                                          fluid );
+    }
+    else
+    {
+      return isothermalCompositionalMultiphaseBaseKernels::
+               PhaseVolumeFractionKernelFactory::
+               createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                          m_numPhases,
+                                                          dataGroup,
+                                                          fluid );
+    }
+  }
 }
 
 void CompositionalMultiphaseBase::updateFluidModel( ObjectManagerBase & dataGroup ) const
@@ -755,22 +839,43 @@ void CompositionalMultiphaseBase::updateCompAmount( ElementSubRegionBase & subRe
 {
   GEOS_MARK_FUNCTION;
 
+  arrayView2d< real64, compflow::USD_COMP > const compAmount = subRegion.getField< fields::flow::compAmount >();
+  arrayView1d< real64 const > const volume = subRegion.getElementVolume();
   string const & solidName = subRegion.template getReference< string >( viewKeyStruct::solidNamesString() );
   CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName );
   arrayView2d< real64 const > const porosity = porousMaterial.getPorosity();
-  arrayView1d< real64 const > const volume = subRegion.getElementVolume();
-  arrayView2d< real64 const, compflow::USD_COMP > const compDens = subRegion.getField< fields::flow::globalCompDensity >();
-  arrayView2d< real64, compflow::USD_COMP > const compAmount = subRegion.getField< fields::flow::compAmount >();
 
   integer const numComp = m_numComponents;
 
-  forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
+  if( m_useZFormulation )
   {
-    for( integer ic = 0; ic < numComp; ++ic )
+    arrayView2d< real64 const, compflow::USD_COMP > const compFrac = subRegion.getField< fields::flow::globalCompFraction >();
+    // access total density stored in the fluid
+    string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() );
+    MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName );
+    arrayView2d< real64 const, multifluid::USD_FLUID > const totalDens = fluid.totalDensity();
+
+    forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
     {
-      compAmount[ei][ic] = porosity[ei][0] * volume[ei] * compDens[ei][ic];
-    }
-  } );
+      for( integer ic = 0; ic < numComp; ++ic )
+      {
+        // m = phi*V*z_c*rho_T
+        compAmount[ei][ic] = porosity[ei][0] * volume[ei] * compFrac[ei][ic] * totalDens[ei][0];
+      }
+    } );
+  }
+  else
+  {
+    arrayView2d< real64 const, compflow::USD_COMP > const compDens = subRegion.getField< fields::flow::globalCompDensity >();
+
+    forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
+    {
+      for( integer ic = 0; ic < numComp; ++ic )
+      {
+        compAmount[ei][ic] = porosity[ei][0] * volume[ei] * compDens[ei][ic];
+      }
+    } );
+  }
 }
 
 void CompositionalMultiphaseBase::updateEnergy( ElementSubRegionBase & subRegion ) const
@@ -824,16 +929,22 @@ real64 CompositionalMultiphaseBase::updateFluidState( ElementSubRegionBase & sub
 {
   GEOS_MARK_FUNCTION;
 
-  updateGlobalComponentFraction( subRegion );
+  real64 maxDeltaPhaseVolFrac;
+
+  if( !m_useZFormulation )
+  {
+    // For p, rho_c as the primary unknowns
+    updateGlobalComponentFraction( subRegion );
+  }
   updateFluidModel( subRegion );
   updateCompAmount( subRegion );
-  real64 const maxDeltaPhaseVolFrac = updatePhaseVolumeFraction( subRegion );
+  maxDeltaPhaseVolFrac = updatePhaseVolumeFraction( subRegion );
   updateRelPermModel( subRegion );
   updatePhaseMobility( subRegion );
   updateCapPressureModel( subRegion );
-
   // note1: for now, thermal conductivity is treated explicitly, so no update here
   // note2: for now, diffusion and dispersion are also treated explicitly
+
   return maxDeltaPhaseVolFrac;
 }
 
@@ -857,26 +968,37 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh,
     string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() );
     MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName );
     arrayView2d< real64 const, multifluid::USD_FLUID > const totalDens = fluid.totalDensity();
-    arrayView2d< real64 const, compflow::USD_COMP > const compFrac =
-      subRegion.getField< fields::flow::globalCompFraction >();
-    arrayView2d< real64, compflow::USD_COMP > const compDens =
-      subRegion.getField< fields::flow::globalCompDensity >();
 
-    forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
+    if( !m_useZFormulation )
     {
-      for( integer ic = 0; ic < numComp; ++ic )
+      arrayView2d< real64 const, compflow::USD_COMP > const compFrac =
+        subRegion.getField< fields::flow::globalCompFraction >();
+      arrayView2d< real64, compflow::USD_COMP > const compDens =
+        subRegion.getField< fields::flow::globalCompDensity >();
+
+      forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
       {
-        compDens[ei][ic] = totalDens[ei][0] * compFrac[ei][ic];
-      }
-    } );
+        for( integer ic = 0; ic < numComp; ++ic )
+        {
+          compDens[ei][ic] = totalDens[ei][0] * compFrac[ei][ic];
+        }
+      } );
 
-    // with initial component densities defined - check if they need to be corrected to avoid zero diags etc
-    if( m_allowCompDensChopping )
-    {
-      chopNegativeDensities( subRegion );
+      // with initial component densities defined - check if they need to be corrected to avoid zero diags etc
+      if( m_allowCompDensChopping )
+      {
+        chopNegativeDensities( subRegion );
+      }
     }
   } );
 
+  // check if comp fractions need to be corrected to avoid zero diags etc
+  if( m_useZFormulation && m_allowCompDensChopping )
+  {
+    DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+    chopNegativeCompFractions( domain );
+  }
+
   // for some reason CUDA does not want the host_device lambda to be defined inside the generic lambda
   // I need the exact type of the subRegion for updateSolidflowProperties to work well.
   mesh.getElemManager().forElementSubRegions< CellElementSubRegion,
@@ -905,9 +1027,9 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh,
     // - This step depends on phaseVolFraction
     RelativePermeabilityBase & relPerm =
       getConstitutiveModel< RelativePermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::relPermNamesString() ) );
-    relPerm.saveConvergedPhaseVolFractionState( phaseVolFrac ); // this needs to happen before calling updateRelPermModel
+    relPerm.saveConvergedPhaseVolFractionState( phaseVolFrac );   // this needs to happen before calling updateRelPermModel
     updateRelPermModel( subRegion );
-    relPerm.saveConvergedState(); // this needs to happen after calling updateRelPermModel
+    relPerm.saveConvergedState();   // this needs to happen after calling updateRelPermModel
 
     string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() );
     MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName );
@@ -937,7 +1059,7 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh,
 
       CapillaryPressureBase const & capPressure =
         getConstitutiveModel< CapillaryPressureBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::capPressureNamesString() ) );
-      capPressure.initializeRockState( porosity, permeability ); // this needs to happen before calling updateCapPressureModel
+      capPressure.initializeRockState( porosity, permeability );   // this needs to happen before calling updateCapPressureModel
       updateCapPressureModel( subRegion );
     }
 
@@ -1069,7 +1191,7 @@ void CompositionalMultiphaseBase::computeHydrostaticEquilibrium( DomainPartition
       real64 const equilTolerance = fs.getEquilibrationTolerance();
       real64 const datumElevation = fs.getDatumElevation();
       real64 const datumPressure = fs.getDatumPressure();
-      string const initPhaseName = fs.getInitPhaseName(); // will go away when GOC/WOC are implemented
+      string const initPhaseName = fs.getInitPhaseName();   // will go away when GOC/WOC are implemented
 
       localIndex const equilIndex = equilNameToEquilId.at( fs.getName() );
       real64 const minElevation = LvArray::math::min( globalMinElevation[equilIndex], datumElevation );
@@ -1077,7 +1199,7 @@ void CompositionalMultiphaseBase::computeHydrostaticEquilibrium( DomainPartition
       real64 const elevationIncrement = LvArray::math::min( fs.getElevationIncrement(), maxElevation - minElevation );
       localIndex const numPointsInTable = ( elevationIncrement > 0 ) ? std::ceil( (maxElevation - minElevation) / elevationIncrement ) + 1 : 1;
 
-      real64 const eps = 0.1 * (maxElevation - minElevation); // we add a small buffer to only log in the pathological cases
+      real64 const eps = 0.1 * (maxElevation - minElevation);   // we add a small buffer to only log in the pathological cases
       GEOS_LOG_RANK_0_IF( ( (datumElevation > globalMaxElevation[equilIndex]+eps)  || (datumElevation < globalMinElevation[equilIndex]-eps) ),
                           getCatalogName() << " " << getDataContext() <<
                           ": By looking at the elevation of the cell centers in this model, GEOS found that " <<
@@ -1116,7 +1238,7 @@ void CompositionalMultiphaseBase::computeHydrostaticEquilibrium( DomainPartition
       auto itRegionFilter = regionFilter.find( region.getName() );
       if( itRegionFilter == regionFilter.end() )
       {
-        return; // the region is not in target, there is nothing to do
+        return;   // the region is not in target, there is nothing to do
       }
       string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() );
       MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName );
@@ -1254,7 +1376,7 @@ void CompositionalMultiphaseBase::initializePostInitialConditionsPreSubGroups()
                                                                arrayView1d< string const > const & regionNames )
   {
     FieldIdentifiers fieldsToBeSync;
-    fieldsToBeSync.addElementFields( { fields::flow::globalCompDensity::key() },
+    fieldsToBeSync.addElementFields( {m_useZFormulation ? fields::flow::globalCompFraction::key() : fields::flow::globalCompDensity::key() },
                                      regionNames );
 
     CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false );
@@ -1319,10 +1441,11 @@ void CompositionalMultiphaseBase::assembleSystem( real64 const GEOS_UNUSED_PARAM
 {
   GEOS_MARK_FUNCTION;
 
-  assembleAccumulationAndVolumeBalanceTerms( domain,
-                                             dofManager,
-                                             localMatrix,
-                                             localRhs );
+  // Accumulation and other local terms
+  assembleLocalTerms( domain,
+                      dofManager,
+                      localMatrix,
+                      localRhs );
 
   if( m_isJumpStabilized )
   {
@@ -1342,10 +1465,10 @@ void CompositionalMultiphaseBase::assembleSystem( real64 const GEOS_UNUSED_PARAM
   }
 }
 
-void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( DomainPartition & domain,
-                                                                             DofManager const & dofManager,
-                                                                             CRSMatrixView< real64, globalIndex const > const & localMatrix,
-                                                                             arrayView1d< real64 > const & localRhs ) const
+void CompositionalMultiphaseBase::assembleLocalTerms( DomainPartition & domain,
+                                                      DofManager const & dofManager,
+                                                      CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                                                      arrayView1d< real64 > const & localRhs ) const
 {
   GEOS_MARK_FUNCTION;
 
@@ -1364,14 +1487,16 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom
       MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName );
       CoupledSolidBase const & solid = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName );
 
-      if( m_isThermal )
+      if( m_useZFormulation )
       {
-        thermalCompositionalMultiphaseBaseKernels::
-          AccumulationKernelFactory::
+        // isothermal for now
+        isothermalCompositionalMultiphaseBaseKernels::
+          AccumulationZFormulationKernelFactory::
           createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
                                                      m_numPhases,
                                                      dofManager.rankOffset(),
                                                      m_useTotalMassEquation,
+                                                     m_useSimpleAccumulation,
                                                      dofKey,
                                                      subRegion,
                                                      fluid,
@@ -1381,19 +1506,37 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom
       }
       else
       {
-        isothermalCompositionalMultiphaseBaseKernels::
-          AccumulationKernelFactory::
-          createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
-                                                     m_numPhases,
-                                                     dofManager.rankOffset(),
-                                                     m_useTotalMassEquation,
-                                                     m_useSimpleAccumulation,
-                                                     dofKey,
-                                                     subRegion,
-                                                     fluid,
-                                                     solid,
-                                                     localMatrix,
-                                                     localRhs );
+        if( m_isThermal )
+        {
+          thermalCompositionalMultiphaseBaseKernels::
+            AccumulationKernelFactory::
+            createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                       m_numPhases,
+                                                       dofManager.rankOffset(),
+                                                       m_useTotalMassEquation,
+                                                       dofKey,
+                                                       subRegion,
+                                                       fluid,
+                                                       solid,
+                                                       localMatrix,
+                                                       localRhs );
+        }
+        else
+        {
+          isothermalCompositionalMultiphaseBaseKernels::
+            AccumulationKernelFactory::
+            createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                       m_numPhases,
+                                                       dofManager.rankOffset(),
+                                                       m_useTotalMassEquation,
+                                                       m_useSimpleAccumulation,
+                                                       dofKey,
+                                                       subRegion,
+                                                       fluid,
+                                                       solid,
+                                                       localMatrix,
+                                                       localRhs );
+        }
       }
     } );
   } );
@@ -1574,7 +1717,7 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time,
           return;
         }
 
-        real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor here!
+        real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor;   // scale the contribution by the sizeScalingFactor here!
         massProd += rhsValue;
         if( useTotalMassEquation > 0 )
         {
@@ -1583,7 +1726,7 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time,
           localRhs[totalMassBalanceRow] += rhsValue;
           if( fluidComponentId < numFluidComponents - 1 )
           {
-            globalIndex const compMassBalanceRow = totalMassBalanceRow + fluidComponentId + 1; // component mass bal equations are shifted
+            globalIndex const compMassBalanceRow = totalMassBalanceRow + fluidComponentId + 1;   // component mass bal equations are shifted
             localRhs[compMassBalanceRow] += rhsValue;
           }
         }
@@ -1711,7 +1854,7 @@ bool CompositionalMultiphaseBase::validateDirichletBC( DomainPartition & domain,
         bcConsistent = false;
         GEOS_WARNING( BCMessage::invalidComponentIndex( comp, fs.getName(),
                                                         fields::flow::globalCompFraction::key() ) );
-        return; // can't check next part with invalid component id
+        return;   // can't check next part with invalid component id
       }
 
       ComponentMask< MAX_NC > & compMask = subRegionSetMap[setName];
@@ -1845,46 +1988,93 @@ void CompositionalMultiphaseBase::applyDirichletBC( real64 const time_n,
         subRegion.getReference< array1d< globalIndex > >( dofKey );
       arrayView1d< real64 const > const pres =
         subRegion.getReference< array1d< real64 > >( fields::flow::pressure::key() );
-      arrayView2d< real64 const, compflow::USD_COMP > const compDens =
-        subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( fields::flow::globalCompDensity::key() );
-      arrayView2d< real64 const, multifluid::USD_FLUID > const totalDens = fluid.totalDensity();
 
-      integer const numComp = m_numComponents;
-      forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
+      // for now to avoid lambda function complain
+
+      if( m_useZFormulation )
       {
-        localIndex const ei = targetSet[a];
-        if( ghostRank[ei] >= 0 )
+        //arrayView2d< real64 const, compflow::USD_COMP > const compFrac =
+        //subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( fields::flow::globalCompFraction::key() );
+
+        integer const numComp = m_numComponents;
+        forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
         {
-          return;
-        }
+          localIndex const ei = targetSet[a];
+          if( ghostRank[ei] >= 0 )
+          {
+            return;
+          }
 
-        globalIndex const dofIndex = dofNumber[ei];
-        localIndex const localRow = dofIndex - rankOffset;
-        real64 rhsValue;
+          globalIndex const dofIndex = dofNumber[ei];
+          localIndex const localRow = dofIndex - rankOffset;
+          real64 rhsValue;
 
-        // 4.1. Apply pressure value to the matrix/rhs
-        FieldSpecificationEqual::SpecifyFieldValue( dofIndex,
-                                                    rankOffset,
-                                                    localMatrix,
-                                                    rhsValue,
-                                                    bcPres[ei],
-                                                    pres[ei] );
-        localRhs[localRow] = rhsValue;
+          // 4.1. Apply pressure value to the matrix/rhs
+          FieldSpecificationEqual::SpecifyFieldValue( dofIndex,
+                                                      rankOffset,
+                                                      localMatrix,
+                                                      rhsValue,
+                                                      bcPres[ei],
+                                                      pres[ei] );
+          localRhs[localRow] = rhsValue;
 
-        // 4.2. For each component, apply target global density value
-        for( integer ic = 0; ic < numComp; ++ic )
+          // 4.2. For each component, apply target global density value
+          for( integer ic = 0; ic < numComp; ++ic )
+          {
+            FieldSpecificationEqual::SpecifyFieldValue( dofIndex + ic + 1,
+                                                        rankOffset,
+                                                        localMatrix,
+                                                        rhsValue,
+                                                        compFrac[ei][ic],
+                                                        compFrac[ei][ic] );
+            localRhs[localRow + ic + 1] = rhsValue;
+          }
+        } );
+      }
+      else
+      {
+        arrayView2d< real64 const, compflow::USD_COMP > const compDens =
+          subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( fields::flow::globalCompDensity::key() );
+        arrayView2d< real64 const, multifluid::USD_FLUID > const totalDens = fluid.totalDensity();
+
+        integer const numComp = m_numComponents;
+        forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
         {
-          FieldSpecificationEqual::SpecifyFieldValue( dofIndex + ic + 1,
+          localIndex const ei = targetSet[a];
+          if( ghostRank[ei] >= 0 )
+          {
+            return;
+          }
+
+          globalIndex const dofIndex = dofNumber[ei];
+          localIndex const localRow = dofIndex - rankOffset;
+          real64 rhsValue;
+
+          // 4.1. Apply pressure value to the matrix/rhs
+          FieldSpecificationEqual::SpecifyFieldValue( dofIndex,
                                                       rankOffset,
                                                       localMatrix,
                                                       rhsValue,
-                                                      totalDens[ei][0] * compFrac[ei][ic],
-                                                      compDens[ei][ic] );
-          localRhs[localRow + ic + 1] = rhsValue;
-        }
-      } );
+                                                      bcPres[ei],
+                                                      pres[ei] );
+          localRhs[localRow] = rhsValue;
+
+          // 4.2. For each component, apply target global component fraction value
+          for( integer ic = 0; ic < numComp; ++ic )
+          {
+            FieldSpecificationEqual::SpecifyFieldValue( dofIndex + ic + 1,
+                                                        rankOffset,
+                                                        localMatrix,
+                                                        rhsValue,
+                                                        totalDens[ei][0] * compFrac[ei][ic],
+                                                        compDens[ei][ic] );
+            localRhs[localRow + ic + 1] = rhsValue;
+          }
+        } );
+      }
     } );
 
+
     // 5. Apply temperature to the system
     if( m_isThermal )
     {
@@ -1980,7 +2170,7 @@ void CompositionalMultiphaseBase::keepVariablesConstantDuringInitStep( real64 co
                                                     rankOffset,
                                                     localMatrix,
                                                     rhsValue,
-                                                    pres[ei], // freeze the current pressure value
+                                                    pres[ei],   // freeze the current pressure value
                                                     pres[ei] );
         localRhs[localRow] = rhsValue;
 
@@ -1991,7 +2181,7 @@ void CompositionalMultiphaseBase::keepVariablesConstantDuringInitStep( real64 co
                                                       rankOffset,
                                                       localMatrix,
                                                       rhsValue,
-                                                      temp[ei], // freeze the current temperature value
+                                                      temp[ei],   // freeze the current temperature value
                                                       temp[ei] );
           localRhs[localRow + numComp + 1] = rhsValue;
         }
@@ -2003,7 +2193,7 @@ void CompositionalMultiphaseBase::keepVariablesConstantDuringInitStep( real64 co
                                                       rankOffset,
                                                       localMatrix,
                                                       rhsValue,
-                                                      compDens[ei][ic], // freeze the current component density values
+                                                      compDens[ei][ic],   // freeze the current component density values
                                                       compDens[ei][ic] );
           localRhs[localRow + ic + 1] = rhsValue;
         }
@@ -2016,8 +2206,6 @@ void CompositionalMultiphaseBase::chopNegativeDensities( DomainPartition & domai
 {
   GEOS_MARK_FUNCTION;
 
-  using namespace isothermalCompositionalMultiphaseBaseKernels;
-
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
                                                                MeshLevel & mesh,
                                                                arrayView1d< string const > const & regionNames )
@@ -2056,24 +2244,14 @@ void CompositionalMultiphaseBase::chopNegativeDensities( ElementSubRegionBase &
   } );
 }
 
-real64 CompositionalMultiphaseBase::setNextDtBasedOnStateChange( real64 const & currentDt,
-                                                                 DomainPartition & domain )
+void CompositionalMultiphaseBase::chopNegativeCompFractions( DomainPartition & domain )
 {
-  if( m_targetRelativePresChange >= 1.0 &&
-      m_targetPhaseVolFracChange >= 1.0 &&
-      m_targetRelativeCompDensChange >= 1.0 &&
-      ( !m_isThermal || m_targetRelativeTempChange >= 1.0 ) )
-  {
-    return LvArray::NumericLimits< real64 >::max;
-  }
+  GEOS_MARK_FUNCTION;
 
-  real64 maxRelativePresChange = 0.0;
-  real64 maxRelativeTempChange = 0.0;
-  real64 maxAbsolutePhaseVolFracChange = 0.0;
-  real64 maxRelativeCompDensChange = 0.0;
+  using namespace isothermalCompositionalMultiphaseBaseKernels;
 
-  integer const numPhase = m_numPhases;
   integer const numComp = m_numComponents;
+  real64 const minCompFrac = m_minCompFrac;
 
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
                                                                MeshLevel & mesh,
@@ -2085,100 +2263,276 @@ real64 CompositionalMultiphaseBase::setNextDtBasedOnStateChange( real64 const &
     {
       arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
 
-      arrayView1d< real64 const > const pres = subRegion.getField< fields::flow::pressure >();
-      arrayView1d< real64 const > const pres_n = subRegion.getField< fields::flow::pressure_n >();
-      arrayView1d< real64 const > const temp = subRegion.getField< fields::flow::temperature >();
-      arrayView1d< real64 const > const temp_n = subRegion.getField< fields::flow::temperature_n >();
-      arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac =
-        subRegion.getField< fields::flow::phaseVolumeFraction >();
-      arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac_n =
-        subRegion.getField< fields::flow::phaseVolumeFraction_n >();
-      arrayView2d< real64 const, compflow::USD_COMP > const compDens =
-        subRegion.getField< fields::flow::globalCompDensity >();
-      arrayView2d< real64, compflow::USD_COMP > const compDens_n =
-        subRegion.getField< fields::flow::globalCompDensity_n >();
-
-      RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxPresChange( 0.0 );
-      RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxTempChange( 0.0 );
-      RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxPhaseVolFracChange( 0.0 );
-      RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxCompDensChange( 0.0 );
+      arrayView2d< real64, compflow::USD_COMP > const compFrac =
+        subRegion.getField< fields::flow::globalCompFraction >();
 
       forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
       {
         if( ghostRank[ei] < 0 )
         {
-          // switch from relative to absolute when values less than 1
-          subRegionMaxPresChange.max( LvArray::math::abs( pres[ei] - pres_n[ei] ) / LvArray::math::max( LvArray::math::abs( pres_n[ei] ), 1.0 ) );
-          subRegionMaxTempChange.max( LvArray::math::abs( temp[ei] - temp_n[ei] ) / LvArray::math::max( LvArray::math::abs( temp_n[ei] ), 1.0 ) );
-          for( integer ip = 0; ip < numPhase; ++ip )
-          {
-            subRegionMaxPhaseVolFracChange.max( LvArray::math::abs( phaseVolFrac[ei][ip] - phaseVolFrac_n[ei][ip] ) );
-          }
           for( integer ic = 0; ic < numComp; ++ic )
           {
-            subRegionMaxCompDensChange.max( LvArray::math::abs( compDens[ei][ic] - compDens_n[ei][ic] ) / LvArray::math::max( LvArray::math::abs( compDens_n[ei][ic] ), 1.0 ) );
+            if( compFrac[ei][ic] < minCompFrac )
+            {
+              compFrac[ei][ic] = minCompFrac;
+            }
           }
         }
       } );
+    } );
+  } );
+}
+
+real64 CompositionalMultiphaseBase::setNextDtBasedOnStateChange( real64 const & currentDt,
+                                                                 DomainPartition & domain )
+{
+  // TODO: put in a seprate function or put flags to avoid duplication of code
+  if( m_useZFormulation )
+  {
+    if( m_targetRelativePresChange >= 1.0 &&
+        m_targetPhaseVolFracChange >= 1.0 &&
+        m_targetCompFracChange >= 1.0 &&
+        ( !m_isThermal || m_targetRelativeTempChange >= 1.0 ) )
+    {
+      return LvArray::NumericLimits< real64 >::max;
+    }
+
+    real64 maxRelativePresChange = 0.0;
+    real64 maxRelativeTempChange = 0.0;
+    real64 maxAbsolutePhaseVolFracChange = 0.0;
+    real64 maxCompFracChange = 0.0;
+
+    integer const numPhase = m_numPhases;
+    integer const numComp = m_numComponents;
+
+    forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+                                                                 MeshLevel & mesh,
+                                                                 arrayView1d< string const > const & regionNames )
+    {
+      mesh.getElemManager().forElementSubRegions( regionNames,
+                                                  [&]( localIndex const,
+                                                       ElementSubRegionBase & subRegion )
+      {
+        arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
+
+        arrayView1d< real64 const > const pres = subRegion.getField< fields::flow::pressure >();
+        arrayView1d< real64 const > const pres_n = subRegion.getField< fields::flow::pressure_n >();
+        arrayView1d< real64 const > const temp = subRegion.getField< fields::flow::temperature >();
+        arrayView1d< real64 const > const temp_n = subRegion.getField< fields::flow::temperature_n >();
+        arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac =
+          subRegion.getField< fields::flow::phaseVolumeFraction >();
+        arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac_n =
+          subRegion.getField< fields::flow::phaseVolumeFraction_n >();
+        arrayView2d< real64 const, compflow::USD_COMP > const compFrac =
+          subRegion.getField< fields::flow::globalCompFraction >();
+        arrayView2d< real64, compflow::USD_COMP > const compFrac_n =
+          subRegion.getField< fields::flow::globalCompFraction_n >();
+
+        RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxPresChange( 0.0 );
+        RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxTempChange( 0.0 );
+        RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxPhaseVolFracChange( 0.0 );
+        RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxCompFracChange( 0.0 );
+
+        forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
+        {
+          if( ghostRank[ei] < 0 )
+          {
+            // switch from relative to absolute when values less than 1
+            subRegionMaxPresChange.max( LvArray::math::abs( pres[ei] - pres_n[ei] ) / LvArray::math::max( LvArray::math::abs( pres_n[ei] ), 1.0 ) );
+            subRegionMaxTempChange.max( LvArray::math::abs( temp[ei] - temp_n[ei] ) / LvArray::math::max( LvArray::math::abs( temp_n[ei] ), 1.0 ) );
+            for( integer ip = 0; ip < numPhase; ++ip )
+            {
+              subRegionMaxPhaseVolFracChange.max( LvArray::math::abs( phaseVolFrac[ei][ip] - phaseVolFrac_n[ei][ip] ) );
+            }
+            for( integer ic = 0; ic < numComp; ++ic )
+            {
+              subRegionMaxCompFracChange.max( LvArray::math::abs( compFrac[ei][ic] - compFrac_n[ei][ic] ) );
+            }
+          }
+        } );
 
-      maxRelativePresChange = LvArray::math::max( maxRelativePresChange, subRegionMaxPresChange.get() );
-      maxRelativeTempChange = LvArray::math::max( maxRelativeTempChange, subRegionMaxTempChange.get() );
-      maxAbsolutePhaseVolFracChange = LvArray::math::max( maxAbsolutePhaseVolFracChange, subRegionMaxPhaseVolFracChange.get() );
-      maxRelativeCompDensChange = LvArray::math::max( maxRelativeCompDensChange, subRegionMaxCompDensChange.get() );
+        maxRelativePresChange = LvArray::math::max( maxRelativePresChange, subRegionMaxPresChange.get() );
+        maxRelativeTempChange = LvArray::math::max( maxRelativeTempChange, subRegionMaxTempChange.get() );
+        maxAbsolutePhaseVolFracChange = LvArray::math::max( maxAbsolutePhaseVolFracChange, subRegionMaxPhaseVolFracChange.get() );
+        maxCompFracChange = LvArray::math::max( maxCompFracChange, subRegionMaxCompFracChange.get() );
 
+      } );
     } );
-  } );
 
-  maxRelativePresChange = MpiWrapper::max( maxRelativePresChange );
-  maxAbsolutePhaseVolFracChange = MpiWrapper::max( maxAbsolutePhaseVolFracChange );
+    maxRelativePresChange = MpiWrapper::max( maxRelativePresChange );
+    maxAbsolutePhaseVolFracChange = MpiWrapper::max( maxAbsolutePhaseVolFracChange );
 
-  GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max relative pressure change during time step = {} %",
-                                                           getName(), GEOS_FMT( "{:.{}f}", 100*maxRelativePresChange, 3 ) ) );
-  GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max absolute phase volume fraction change during time step = {}",
-                                                           getName(), GEOS_FMT( "{:.{}f}", maxAbsolutePhaseVolFracChange, 3 ) ) );
+    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max relative pressure change during time step = {} %",
+                                                             getName(), GEOS_FMT( "{:.{}f}", 100*maxRelativePresChange, 3 ) ) );
+    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max absolute phase volume fraction change during time step = {}",
+                                                             getName(), GEOS_FMT( "{:.{}f}", maxAbsolutePhaseVolFracChange, 3 ) ) );
 
-  if( m_targetRelativeCompDensChange < LvArray::NumericLimits< real64 >::max )
-  {
-    maxRelativeCompDensChange = MpiWrapper::max( maxRelativeCompDensChange );
-    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max relative component density change during time step = {} %",
-                                                             getName(), GEOS_FMT( "{:.{}f}", 100*maxRelativeCompDensChange, 3 ) ) );
-  }
+    if( m_targetCompFracChange < LvArray::NumericLimits< real64 >::max )
+    {
+      maxCompFracChange = MpiWrapper::max( maxCompFracChange );
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max component fraction change during time step = {} %",
+                                                               getName(), GEOS_FMT( "{:.{}f}", 100*maxCompFracChange, 3 ) ) );
+    }
 
-  if( m_isThermal )
-  {
-    maxRelativeTempChange = MpiWrapper::max( maxRelativeTempChange );
-    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max relative temperature change during time step = {} %",
-                                                             getName(), GEOS_FMT( "{:.{}f}", 100*maxRelativeTempChange, 3 ) ) );
-  }
+    if( m_isThermal )
+    {
+      maxRelativeTempChange = MpiWrapper::max( maxRelativeTempChange );
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max relative temperature change during time step = {} %",
+                                                               getName(), GEOS_FMT( "{:.{}f}", 100*maxRelativeTempChange, 3 ) ) );
+    }
 
-  real64 const eps = LvArray::NumericLimits< real64 >::epsilon;
-
-  real64 const nextDtPressure = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetRelativePresChange
-                                / std::max( eps, maxRelativePresChange + m_solutionChangeScalingFactor * m_targetRelativePresChange );
-  if( m_nonlinearSolverParameters.getLogLevel() > 0 )
-    GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on pressure change = {}", getName(), nextDtPressure ));
-  real64 const nextDtPhaseVolFrac = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetPhaseVolFracChange
-                                    / std::max( eps, maxAbsolutePhaseVolFracChange + m_solutionChangeScalingFactor * m_targetPhaseVolFracChange );
-  if( m_nonlinearSolverParameters.getLogLevel() > 0 )
-    GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on phase volume fraction change = {}", getName(), nextDtPhaseVolFrac ));
-  real64 nextDtCompDens = LvArray::NumericLimits< real64 >::max;
-  if( m_targetRelativeCompDensChange < LvArray::NumericLimits< real64 >::max )
-  {
-    nextDtCompDens = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetRelativeCompDensChange
-                     / std::max( eps, maxRelativeCompDensChange + m_solutionChangeScalingFactor * m_targetRelativeCompDensChange );
+    real64 const eps = LvArray::NumericLimits< real64 >::epsilon;
+
+    real64 const nextDtPressure = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetRelativePresChange
+                                  / std::max( eps, maxRelativePresChange + m_solutionChangeScalingFactor * m_targetRelativePresChange );
+    if( m_nonlinearSolverParameters.getLogLevel() > 0 )
+      GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on pressure change = {}", getName(), nextDtPressure ));
+    real64 const nextDtPhaseVolFrac = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetPhaseVolFracChange
+                                      / std::max( eps, maxAbsolutePhaseVolFracChange + m_solutionChangeScalingFactor * m_targetPhaseVolFracChange );
     if( m_nonlinearSolverParameters.getLogLevel() > 0 )
-      GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on component density change = {}", getName(), nextDtCompDens ));
+      GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on phase volume fraction change = {}", getName(), nextDtPhaseVolFrac ));
+    real64 nextDtCompFrac = LvArray::NumericLimits< real64 >::max;
+    if( m_targetCompFracChange < LvArray::NumericLimits< real64 >::max )
+    {
+      nextDtCompFrac = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetCompFracChange
+                       / std::max( eps, maxCompFracChange + m_solutionChangeScalingFactor * m_targetCompFracChange );
+      if( m_nonlinearSolverParameters.getLogLevel() > 0 )
+        GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on component fraction change = {}", getName(), nextDtCompFrac ));
+    }
+    real64 nextDtTemperature = LvArray::NumericLimits< real64 >::max;
+    if( m_isThermal )
+    {
+      nextDtTemperature = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetRelativeTempChange
+                          / std::max( eps, maxRelativeTempChange + m_solutionChangeScalingFactor * m_targetRelativeTempChange );
+      if( m_nonlinearSolverParameters.getLogLevel() > 0 )
+        GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on temperature change = {}", getName(), nextDtPhaseVolFrac ));
+    }
+
+    return std::min( std::min( nextDtPressure, std::min( nextDtPhaseVolFrac, nextDtCompFrac ) ), nextDtTemperature );
   }
-  real64 nextDtTemperature = LvArray::NumericLimits< real64 >::max;
-  if( m_isThermal )
+  else
   {
-    nextDtTemperature = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetRelativeTempChange
-                        / std::max( eps, maxRelativeTempChange + m_solutionChangeScalingFactor * m_targetRelativeTempChange );
+    if( m_targetRelativePresChange >= 1.0 &&
+        m_targetPhaseVolFracChange >= 1.0 &&
+        m_targetRelativeCompDensChange >= 1.0 &&
+        ( !m_isThermal || m_targetRelativeTempChange >= 1.0 ) )
+    {
+      return LvArray::NumericLimits< real64 >::max;
+    }
+
+    real64 maxRelativePresChange = 0.0;
+    real64 maxRelativeTempChange = 0.0;
+    real64 maxAbsolutePhaseVolFracChange = 0.0;
+    real64 maxRelativeCompDensChange = 0.0;
+
+    integer const numPhase = m_numPhases;
+    integer const numComp = m_numComponents;
+
+    forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+                                                                 MeshLevel & mesh,
+                                                                 arrayView1d< string const > const & regionNames )
+    {
+      mesh.getElemManager().forElementSubRegions( regionNames,
+                                                  [&]( localIndex const,
+                                                       ElementSubRegionBase & subRegion )
+      {
+        arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
+
+        arrayView1d< real64 const > const pres = subRegion.getField< fields::flow::pressure >();
+        arrayView1d< real64 const > const pres_n = subRegion.getField< fields::flow::pressure_n >();
+        arrayView1d< real64 const > const temp = subRegion.getField< fields::flow::temperature >();
+        arrayView1d< real64 const > const temp_n = subRegion.getField< fields::flow::temperature_n >();
+        arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac =
+          subRegion.getField< fields::flow::phaseVolumeFraction >();
+        arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac_n =
+          subRegion.getField< fields::flow::phaseVolumeFraction_n >();
+        arrayView2d< real64 const, compflow::USD_COMP > const compDens =
+          subRegion.getField< fields::flow::globalCompDensity >();
+        arrayView2d< real64, compflow::USD_COMP > const compDens_n =
+          subRegion.getField< fields::flow::globalCompDensity_n >();
+
+        RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxPresChange( 0.0 );
+        RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxTempChange( 0.0 );
+        RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxPhaseVolFracChange( 0.0 );
+        RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxCompDensChange( 0.0 );
+
+        forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
+        {
+          if( ghostRank[ei] < 0 )
+          {
+            // switch from relative to absolute when values less than 1
+            subRegionMaxPresChange.max( LvArray::math::abs( pres[ei] - pres_n[ei] ) / LvArray::math::max( LvArray::math::abs( pres_n[ei] ), 1.0 ) );
+            subRegionMaxTempChange.max( LvArray::math::abs( temp[ei] - temp_n[ei] ) / LvArray::math::max( LvArray::math::abs( temp_n[ei] ), 1.0 ) );
+            for( integer ip = 0; ip < numPhase; ++ip )
+            {
+              subRegionMaxPhaseVolFracChange.max( LvArray::math::abs( phaseVolFrac[ei][ip] - phaseVolFrac_n[ei][ip] ) );
+            }
+            for( integer ic = 0; ic < numComp; ++ic )
+            {
+              subRegionMaxCompDensChange.max( LvArray::math::abs( compDens[ei][ic] - compDens_n[ei][ic] ) / LvArray::math::max( LvArray::math::abs( compDens_n[ei][ic] ), 1.0 ) );
+            }
+          }
+        } );
+
+        maxRelativePresChange = LvArray::math::max( maxRelativePresChange, subRegionMaxPresChange.get() );
+        maxRelativeTempChange = LvArray::math::max( maxRelativeTempChange, subRegionMaxTempChange.get() );
+        maxAbsolutePhaseVolFracChange = LvArray::math::max( maxAbsolutePhaseVolFracChange, subRegionMaxPhaseVolFracChange.get() );
+        maxRelativeCompDensChange = LvArray::math::max( maxRelativeCompDensChange, subRegionMaxCompDensChange.get() );
+
+      } );
+    } );
+
+    maxRelativePresChange = MpiWrapper::max( maxRelativePresChange );
+    maxAbsolutePhaseVolFracChange = MpiWrapper::max( maxAbsolutePhaseVolFracChange );
+
+    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max relative pressure change during time step = {} %",
+                                                             getName(), GEOS_FMT( "{:.{}f}", 100*maxRelativePresChange, 3 ) ) );
+    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max absolute phase volume fraction change during time step = {}",
+                                                             getName(), GEOS_FMT( "{:.{}f}", maxAbsolutePhaseVolFracChange, 3 ) ) );
+
+    if( m_targetRelativeCompDensChange < LvArray::NumericLimits< real64 >::max )
+    {
+      maxRelativeCompDensChange = MpiWrapper::max( maxRelativeCompDensChange );
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max relative component density change during time step = {} %",
+                                                               getName(), GEOS_FMT( "{:.{}f}", 100*maxRelativeCompDensChange, 3 ) ) );
+    }
+
+    if( m_isThermal )
+    {
+      maxRelativeTempChange = MpiWrapper::max( maxRelativeTempChange );
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: max relative temperature change during time step = {} %",
+                                                               getName(), GEOS_FMT( "{:.{}f}", 100*maxRelativeTempChange, 3 ) ) );
+    }
+
+    real64 const eps = LvArray::NumericLimits< real64 >::epsilon;
+
+    real64 const nextDtPressure = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetRelativePresChange
+                                  / std::max( eps, maxRelativePresChange + m_solutionChangeScalingFactor * m_targetRelativePresChange );
     if( m_nonlinearSolverParameters.getLogLevel() > 0 )
-      GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on temperature change = {}", getName(), nextDtPhaseVolFrac ));
-  }
+      GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on pressure change = {}", getName(), nextDtPressure ));
+    real64 const nextDtPhaseVolFrac = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetPhaseVolFracChange
+                                      / std::max( eps, maxAbsolutePhaseVolFracChange + m_solutionChangeScalingFactor * m_targetPhaseVolFracChange );
+    if( m_nonlinearSolverParameters.getLogLevel() > 0 )
+      GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on phase volume fraction change = {}", getName(), nextDtPhaseVolFrac ));
+    real64 nextDtCompDens = LvArray::NumericLimits< real64 >::max;
+    if( m_targetRelativeCompDensChange < LvArray::NumericLimits< real64 >::max )
+    {
+      nextDtCompDens = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetRelativeCompDensChange
+                       / std::max( eps, maxRelativeCompDensChange + m_solutionChangeScalingFactor * m_targetRelativeCompDensChange );
+      if( m_nonlinearSolverParameters.getLogLevel() > 0 )
+        GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on component density change = {}", getName(), nextDtCompDens ));
+    }
+    real64 nextDtTemperature = LvArray::NumericLimits< real64 >::max;
+    if( m_isThermal )
+    {
+      nextDtTemperature = currentDt * ( 1.0 + m_solutionChangeScalingFactor ) * m_targetRelativeTempChange
+                          / std::max( eps, maxRelativeTempChange + m_solutionChangeScalingFactor * m_targetRelativeTempChange );
+      if( m_nonlinearSolverParameters.getLogLevel() > 0 )
+        GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on temperature change = {}", getName(), nextDtPhaseVolFrac ));
+    }
 
-  return std::min( std::min( nextDtPressure, std::min( nextDtPhaseVolFrac, nextDtCompDens ) ), nextDtTemperature );
+    return std::min( std::min( nextDtPressure, std::min( nextDtPhaseVolFrac, nextDtCompDens ) ), nextDtTemperature );
+  }
 }
 
 real64 CompositionalMultiphaseBase::setNextDtBasedOnCFL( const geos::real64 & currentDt, geos::DomainPartition & domain )
@@ -2366,11 +2720,23 @@ void CompositionalMultiphaseBase::resetStateToBeginningOfStep( DomainPartition &
         subRegion.template getField< fields::flow::pressure_n >();
       pres.setValues< parallelDevicePolicy<> >( pres_n );
 
-      arrayView2d< real64, compflow::USD_COMP > const & compDens =
-        subRegion.template getField< fields::flow::globalCompDensity >();
-      arrayView2d< real64 const, compflow::USD_COMP > const & compDens_n =
-        subRegion.template getField< fields::flow::globalCompDensity_n >();
-      compDens.setValues< parallelDevicePolicy<> >( compDens_n );
+      if( m_useZFormulation )
+      {
+        arrayView2d< real64, compflow::USD_COMP > const & compFrac =
+          subRegion.template getField< fields::flow::globalCompFraction >();
+        arrayView2d< real64 const, compflow::USD_COMP > const & compFrac_n =
+          subRegion.template getField< fields::flow::globalCompFraction_n >();
+        compFrac.setValues< parallelDevicePolicy<> >( compFrac_n );
+      }
+      else
+      {
+        arrayView2d< real64, compflow::USD_COMP > const & compDens =
+          subRegion.template getField< fields::flow::globalCompDensity >();
+        arrayView2d< real64 const, compflow::USD_COMP > const & compDens_n =
+          subRegion.template getField< fields::flow::globalCompDensity_n >();
+        compDens.setValues< parallelDevicePolicy<> >( compDens_n );
+      }
+
 
       if( m_isThermal )
       {
@@ -2429,11 +2795,11 @@ void CompositionalMultiphaseBase::implicitStepComplete( real64 const & time,
       CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName );
       if( m_keepVariablesConstantDuringInitStep )
       {
-        porousMaterial.ignoreConvergedState(); // newPorosity <- porosity_n
+        porousMaterial.ignoreConvergedState();   // newPorosity <- porosity_n
       }
       else
       {
-        porousMaterial.saveConvergedState(); // porosity_n <- porosity
+        porousMaterial.saveConvergedState();   // porosity_n <- porosity
       }
 
       // Step 4: save converged state for the relperm model to handle hysteresis
@@ -2497,24 +2863,41 @@ void CompositionalMultiphaseBase::saveConvergedState( ElementSubRegionBase & sub
 {
   FlowSolverBase::saveConvergedState( subRegion );
 
-  arrayView2d< real64 const, compflow::USD_COMP > const & compDens =
-    subRegion.template getField< fields::flow::globalCompDensity >();
-  arrayView2d< real64, compflow::USD_COMP > const & compDens_n =
-    subRegion.template getField< fields::flow::globalCompDensity_n >();
-  compDens_n.setValues< parallelDevicePolicy<> >( compDens );
+  if( m_useZFormulation )
+  {
+    arrayView2d< real64 const, compflow::USD_COMP > const & compFrac =
+      subRegion.template getField< fields::flow::globalCompFraction >();
+    arrayView2d< real64, compflow::USD_COMP > const & compFrac_n =
+      subRegion.template getField< fields::flow::globalCompFraction_n >();
+    compFrac_n.setValues< parallelDevicePolicy<> >( compFrac );
+    // may be useful later for sequential poromechanics implementation
+    //if( m_isFixedStressPoromechanicsUpdate )
+    //{
+    //  arrayView2d< real64, compflow::USD_COMP > const & compFrac_k =
+    //    subRegion.template getField< fields::flow::globalCompFraction_k >();
+    //  compFrac_k.setValues< parallelDevicePolicy<> >( compFrac );
+    //}
+  }
+  else
+  {
+    arrayView2d< real64 const, compflow::USD_COMP > const & compDens =
+      subRegion.template getField< fields::flow::globalCompDensity >();
+    arrayView2d< real64, compflow::USD_COMP > const & compDens_n =
+      subRegion.template getField< fields::flow::globalCompDensity_n >();
+    compDens_n.setValues< parallelDevicePolicy<> >( compDens );
+    if( m_isFixedStressPoromechanicsUpdate )
+    {
+      arrayView2d< real64, compflow::USD_COMP > const & compDens_k =
+        subRegion.template getField< fields::flow::globalCompDensity_k >();
+      compDens_k.setValues< parallelDevicePolicy<> >( compDens );
+    }
+  }
 
   arrayView2d< real64 const, compflow::USD_COMP > const & compAmount =
     subRegion.template getField< fields::flow::compAmount >();
   arrayView2d< real64, compflow::USD_COMP > const & compAmount_n =
     subRegion.template getField< fields::flow::compAmount_n >();
   compAmount_n.setValues< parallelDevicePolicy<> >( compAmount );
-
-  if( m_isFixedStressPoromechanicsUpdate )
-  {
-    arrayView2d< real64, compflow::USD_COMP > const & compDens_k =
-      subRegion.template getField< fields::flow::globalCompDensity_k >();
-    compDens_k.setValues< parallelDevicePolicy<> >( compDens );
-  }
 }
 
 void CompositionalMultiphaseBase::saveSequentialIterationState( DomainPartition & domain )
@@ -2559,7 +2942,7 @@ void CompositionalMultiphaseBase::saveSequentialIterationState( DomainPartition
     } );
   } );
 
-  m_sequentialCompDensChange = MpiWrapper::max( maxCompDensChange ); // store to be later used for convergence check
+  m_sequentialCompDensChange = MpiWrapper::max( maxCompDensChange );   // store to be later used for convergence check
 }
 
 void CompositionalMultiphaseBase::updateState( DomainPartition & domain )
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp
index 1cec9ca685..174e2c856d 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp
@@ -191,7 +191,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase
   { return m_useMass ? units::Unit::Mass : units::Unit::Mole; }
 
   /**
-   * @brief assembles the accumulation and volume balance terms for all cells
+   * @brief assembles the accumulation other local terms for all cells
    * @param time_n previous time value
    * @param dt time step
    * @param domain the physical domain object
@@ -199,10 +199,10 @@ class CompositionalMultiphaseBase : public FlowSolverBase
    * @param localMatrix the system matrix
    * @param localRhs the system right-hand side vector
    */
-  void assembleAccumulationAndVolumeBalanceTerms( DomainPartition & domain,
-                                                  DofManager const & dofManager,
-                                                  CRSMatrixView< real64, globalIndex const > const & localMatrix,
-                                                  arrayView1d< real64 > const & localRhs ) const;
+  void assembleLocalTerms( DomainPartition & domain,
+                           DofManager const & dofManager,
+                           CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                           arrayView1d< real64 > const & localRhs ) const;
 
   /**
    * @brief assembles the flux terms for all cells
@@ -234,8 +234,10 @@ class CompositionalMultiphaseBase : public FlowSolverBase
                                DofManager const & dofManager,
                                CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                arrayView1d< real64 > const & localRhs ) const = 0;
+
   /**@}*/
 
+
   struct viewKeyStruct : FlowSolverBase::viewKeyStruct
   {
     static constexpr char const * elemDofFieldString() { return "compositionalVariables"; }
@@ -243,6 +245,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase
     // inputs
 
     static constexpr char const * useMassFlagString() { return "useMass"; }
+    static constexpr char const * useZFormulationFlagString() { return "useZFormulation"; }
     static constexpr char const * relPermNamesString() { return "relPermNames"; }
     static constexpr char const * capPressureNamesString() { return "capPressureNames"; }
     static constexpr char const * diffusionNamesString() { return "diffusionNames"; }
@@ -256,6 +259,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase
     static constexpr char const * targetRelativeTempChangeString() { return "targetRelativeTemperatureChangeInTimeStep"; }
     static constexpr char const * targetPhaseVolFracChangeString() { return "targetPhaseVolFractionChangeInTimeStep"; }
     static constexpr char const * targetRelativeCompDensChangeString() { return "targetRelativeCompDensChangeInTimeStep"; }
+    static constexpr char const * targetCompFracChangeString() { return "targetCompFracChangeInTimeStep"; }
     static constexpr char const * targetFlowCFLString() { return "targetFlowCFL"; }
 
 
@@ -270,6 +274,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase
     static constexpr char const * useSimpleAccumulationString() { return "useSimpleAccumulation"; }
     static constexpr char const * useNewGravityString() { return "useNewGravity"; }
     static constexpr char const * minCompDensString() { return "minCompDens"; }
+    static constexpr char const * minCompFracString() { return "minCompFrac"; }
     static constexpr char const * maxSequentialCompDensChangeString() { return "maxSequentialCompDensChange"; }
     static constexpr char const * minScalingFactorString() { return "minScalingFactor"; }
 
@@ -364,9 +369,14 @@ class CompositionalMultiphaseBase : public FlowSolverBase
    * @param domain the physical domain object
    */
   void chopNegativeDensities( DomainPartition & domain );
-
   void chopNegativeDensities( ElementSubRegionBase & subRegion );
 
+  /**
+   * @brief Sets all the negative component fractions (if any) to zero.
+   * @param domain the physical domain object
+   */
+  void chopNegativeCompFractions( DomainPartition & domain );
+
   virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt,
                                               DomainPartition & domain ) override;
 
@@ -439,6 +449,9 @@ class CompositionalMultiphaseBase : public FlowSolverBase
   /// flag indicating whether mass or molar formulation should be used
   integer m_useMass;
 
+  /// flag indicating whether overall composition (Z) formulation should be used
+  integer m_useZFormulation;
+
   /// flag to determine whether or not to apply capillary pressure
   integer m_hasCapPressure;
 
@@ -475,6 +488,9 @@ class CompositionalMultiphaseBase : public FlowSolverBase
   /// target (relative) change in component density in a time step
   real64 m_targetRelativeCompDensChange;
 
+  /// target (absolute) change in component fraction in a time step
+  real64 m_targetCompFracChange;
+
   /// minimum value of the scaling factor obtained by enforcing maxCompFracChange
   real64 m_minScalingFactor;
 
@@ -493,6 +509,9 @@ class CompositionalMultiphaseBase : public FlowSolverBase
   /// minimum allowed global component density
   real64 m_minCompDens;
 
+  /// minimum allowed global component fraction
+  real64 m_minCompFrac;
+
   /// name of the fluid constitutive model used as a reference for component/phase description
   string m_referenceFluidModelName;
 
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp
index 68faa9d681..3f520747ae 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp
@@ -72,6 +72,23 @@ DECLARE_FIELD( globalCompFraction,
                WRITE_AND_READ,
                "Global component fraction" );
 
+DECLARE_FIELD( globalCompFraction_n,
+               "globalCompFraction_n",
+               array2dLayoutComp,
+               0,
+               NOPLOT,
+               NO_WRITE,
+               "Global component fraction at the previous converged time step" );
+
+// may be needed later for sequential poromechanics implementation
+//DECLARE_FIELD( globalCompFraction_k,
+//               "globalCompFraction_k",
+//               array2dLayoutComp,
+//               0,
+//               NOPLOT,
+//               NO_WRITE,
+//               "Global component fraction updates at the previous sequential iteration" );
+
 DECLARE_FIELD( faceGlobalCompFraction,
                "faceGlobalCompFraction",
                array2dLayoutComp,
@@ -169,6 +186,14 @@ DECLARE_FIELD( globalCompDensityScalingFactor,
                NO_WRITE,
                "Scaling factors for global component densities" );
 
+DECLARE_FIELD( globalCompFractionScalingFactor,
+               "globalCompFractionScalingFactor",
+               array1d< real64 >,
+               1,
+               NOPLOT,
+               NO_WRITE,
+               "Scaling factors for global component fractions" );
+
 DECLARE_FIELD( compAmount,
                "compAmount",
                array2dLayoutComp,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
index 9a64efef23..52cc925d7d 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
@@ -53,6 +53,11 @@
 #include "physicsSolvers/fluidFlow/kernels/compositional/PhaseMobilityKernel.hpp"
 #include "physicsSolvers/fluidFlow/kernels/compositional/ThermalPhaseMobilityKernel.hpp"
 #include "physicsSolvers/fluidFlow/kernels/compositional/AquiferBCKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/zFormulation/FluxComputeZFormulationKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/zFormulation/AccumulationZFormulationKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseMobilityZFormulationKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingZFormulationKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/zFormulation/DirichletFluxComputeZFormulationKernel.hpp"
 
 namespace geos
 {
@@ -106,9 +111,18 @@ void CompositionalMultiphaseFVM::postInputInitialization()
 {
   CompositionalMultiphaseBase::postInputInitialization();
 
-  if( m_scalingType == ScalingType::Local && m_nonlinearSolverParameters.m_lineSearchAction != NonlinearSolverParameters::LineSearchAction::None )
+  if( m_scalingType == ScalingType::Local &&
+      m_nonlinearSolverParameters.m_lineSearchAction != NonlinearSolverParameters::LineSearchAction::None )
   {
-    GEOS_ERROR( GEOS_FMT( "{}: line search is not supported for {} = {}", getName(), viewKeyStruct::scalingTypeString(), EnumStrings< ScalingType >::toString( ScalingType::Local )) );
+    GEOS_ERROR( GEOS_FMT( "{}: line search is not supported for {} = {}",
+                          getName(), viewKeyStruct::scalingTypeString(),
+                          EnumStrings< ScalingType >::toString( ScalingType::Local )) );
+  }
+
+  if( m_useZFormulation && m_dbcParams.useDBC ) // useZFormulation is not compatible with DBC
+  {
+    GEOS_ERROR( GEOS_FMT( "{}: '{}' is not compatible with {}",
+                          getDataContext(), viewKeyStruct::useZFormulationFlagString(), viewKeyStruct::useDBCString() ) );
   }
 }
 
@@ -179,16 +193,20 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt,
       typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper();
 
       // Convective flux
-      if( m_isThermal )
+
+      if( m_useZFormulation )
       {
-        thermalCompositionalMultiphaseFVMKernels::
-          FluxComputeKernelFactory::
+        // isothermal only for now
+        isothermalCompositionalMultiphaseFVMKernels::
+          FluxComputeZFormulationKernelFactory::
           createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
                                                      m_numPhases,
                                                      dofManager.rankOffset(),
                                                      elemDofKey,
                                                      m_hasCapPressure,
                                                      m_useTotalMassEquation,
+                                                     m_useNewGravity,
+                                                     fluxApprox.upwindingParams(),
                                                      getName(),
                                                      mesh.getElemManager(),
                                                      stencilWrapper,
@@ -198,9 +216,9 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt,
       }
       else
       {
-        if( m_dbcParams.useDBC )
+        if( m_isThermal )
         {
-          dissipationCompositionalMultiphaseFVMKernels::
+          thermalCompositionalMultiphaseFVMKernels::
             FluxComputeKernelFactory::
             createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
                                                        m_numPhases,
@@ -213,74 +231,95 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt,
                                                        stencilWrapper,
                                                        dt,
                                                        localMatrix.toViewConstSizes(),
-                                                       localRhs.toView(),
-                                                       m_dbcParams.omega,
-                                                       getNonlinearSolverParameters().m_numNewtonIterations,
-                                                       m_dbcParams.continuation,
-                                                       m_dbcParams.miscible,
-                                                       m_dbcParams.kappamin,
-                                                       m_dbcParams.contMultiplier );
+                                                       localRhs.toView() );
         }
         else
         {
-          isothermalCompositionalMultiphaseFVMKernels::
-            FluxComputeKernelFactory::
-            createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
-                                                       m_numPhases,
-                                                       dofManager.rankOffset(),
-                                                       elemDofKey,
-                                                       m_hasCapPressure,
-                                                       m_useTotalMassEquation,
-                                                       m_useNewGravity,
-                                                       fluxApprox.upwindingParams(),
-                                                       getName(),
-                                                       mesh.getElemManager(),
-                                                       stencilWrapper,
-                                                       dt,
-                                                       localMatrix.toViewConstSizes(),
-                                                       localRhs.toView() );
+          if( m_dbcParams.useDBC )
+          {
+            dissipationCompositionalMultiphaseFVMKernels::
+              FluxComputeKernelFactory::
+              createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                         m_numPhases,
+                                                         dofManager.rankOffset(),
+                                                         elemDofKey,
+                                                         m_hasCapPressure,
+                                                         m_useTotalMassEquation,
+                                                         getName(),
+                                                         mesh.getElemManager(),
+                                                         stencilWrapper,
+                                                         dt,
+                                                         localMatrix.toViewConstSizes(),
+                                                         localRhs.toView(),
+                                                         m_dbcParams.omega,
+                                                         getNonlinearSolverParameters().m_numNewtonIterations,
+                                                         m_dbcParams.continuation,
+                                                         m_dbcParams.miscible,
+                                                         m_dbcParams.kappamin,
+                                                         m_dbcParams.contMultiplier );
+          }
+          else
+          {
+            isothermalCompositionalMultiphaseFVMKernels::
+              FluxComputeKernelFactory::
+              createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                         m_numPhases,
+                                                         dofManager.rankOffset(),
+                                                         elemDofKey,
+                                                         m_hasCapPressure,
+                                                         m_useTotalMassEquation,
+                                                         m_useNewGravity,
+                                                         fluxApprox.upwindingParams(),
+                                                         getName(),
+                                                         mesh.getElemManager(),
+                                                         stencilWrapper,
+                                                         dt,
+                                                         localMatrix.toViewConstSizes(),
+                                                         localRhs.toView() );
+          }
         }
-      }
 
-      // Diffusive and dispersive flux
-      if( m_hasDiffusion || m_hasDispersion )
-      {
+        // Diffusive and dispersive flux
 
-        if( m_isThermal )
+        if( m_hasDiffusion || m_hasDispersion )
         {
-          thermalCompositionalMultiphaseFVMKernels::
-            DiffusionDispersionFluxComputeKernelFactory::
-            createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
-                                                       m_numPhases,
-                                                       dofManager.rankOffset(),
-                                                       elemDofKey,
-                                                       m_hasDiffusion,
-                                                       m_hasDispersion,
-                                                       m_useTotalMassEquation,
-                                                       getName(),
-                                                       mesh.getElemManager(),
-                                                       stencilWrapper,
-                                                       dt,
-                                                       localMatrix.toViewConstSizes(),
-                                                       localRhs.toView() );
-        }
-        else
-        {
-          isothermalCompositionalMultiphaseFVMKernels::
-            DiffusionDispersionFluxComputeKernelFactory::
-            createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
-                                                       m_numPhases,
-                                                       dofManager.rankOffset(),
-                                                       elemDofKey,
-                                                       m_hasDiffusion,
-                                                       m_hasDispersion,
-                                                       m_useTotalMassEquation,
-                                                       getName(),
-                                                       mesh.getElemManager(),
-                                                       stencilWrapper,
-                                                       dt,
-                                                       localMatrix.toViewConstSizes(),
-                                                       localRhs.toView() );
+
+          if( m_isThermal )
+          {
+            thermalCompositionalMultiphaseFVMKernels::
+              DiffusionDispersionFluxComputeKernelFactory::
+              createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                         m_numPhases,
+                                                         dofManager.rankOffset(),
+                                                         elemDofKey,
+                                                         m_hasDiffusion,
+                                                         m_hasDispersion,
+                                                         m_useTotalMassEquation,
+                                                         getName(),
+                                                         mesh.getElemManager(),
+                                                         stencilWrapper,
+                                                         dt,
+                                                         localMatrix.toViewConstSizes(),
+                                                         localRhs.toView() );
+          }
+          else
+          {
+            isothermalCompositionalMultiphaseFVMKernels::
+              DiffusionDispersionFluxComputeKernelFactory::
+              createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                         m_numPhases,
+                                                         dofManager.rankOffset(),
+                                                         elemDofKey,
+                                                         m_hasDiffusion,
+                                                         m_hasDispersion,
+                                                         m_useTotalMassEquation,
+                                                         getName(),
+                                                         mesh.getElemManager(),
+                                                         stencilWrapper,
+                                                         dt,
+                                                         localMatrix.toViewConstSizes(),
+                                                         localRhs.toView() );
+          }
         }
       }
 
@@ -482,14 +521,161 @@ real64 CompositionalMultiphaseFVM::scalingForSystemSolution( DomainPartition & d
 {
   GEOS_MARK_FUNCTION;
 
-  string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
-  real64 scalingFactor = 1.0;
-  real64 minPresScalingFactor = 1.0, minCompDensScalingFactor = 1.0, minTempScalingFactor = 1.0;
+  if( m_useZFormulation )
+  {
+    return scalingForSystemSolutionZFormulation( domain, dofManager, localSolution );
+  }
+  else
+  {
+    string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
+    real64 scalingFactor = 1.0;
+    real64 minPresScalingFactor = 1.0, minCompDensScalingFactor = 1.0, minTempScalingFactor = 1.0;
+
+
+    std::vector< valueAndLocationType > regionDeltaPresMaxLoc;
+    std::vector< valueAndLocationType > regionDeltaTempMaxLoc;
+    std::vector< valueAndLocationType > regionDeltaCompDensMaxLoc;
+
+    forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+                                                                 MeshLevel & mesh,
+                                                                 arrayView1d< string const > const & regionNames )
+    {
+      mesh.getElemManager().forElementSubRegions( regionNames,
+                                                  [&]( localIndex const,
+                                                       ElementSubRegionBase & subRegion )
+      {
+        arrayView1d< globalIndex const > const localToGlobalMap = subRegion.localToGlobalMap();
+        arrayView1d< real64 const > const pressure = subRegion.getField< fields::flow::pressure >();
+        arrayView1d< real64 const > const temperature = subRegion.getField< fields::flow::temperature >();
+        arrayView2d< real64 const, compflow::USD_COMP > const compDens = subRegion.getField< fields::flow::globalCompDensity >();
+        arrayView1d< real64 > pressureScalingFactor = subRegion.getField< fields::flow::pressureScalingFactor >();
+        arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< fields::flow::temperatureScalingFactor >();
+        arrayView1d< real64 > compDensScalingFactor = subRegion.getField< fields::flow::globalCompDensityScalingFactor >();
+
+        const integer temperatureOffset = m_numComponents+1;
+
+        auto const subRegionData =
+          m_isThermal
+    ? thermalCompositionalMultiphaseBaseKernels::
+            SolutionScalingKernelFactory::
+            createAndLaunch< parallelDevicePolicy<> >( m_maxRelativePresChange,
+                                                       m_maxAbsolutePresChange,
+                                                       m_maxRelativeTempChange,
+                                                       m_maxCompFracChange,
+                                                       m_maxRelativeCompDensChange,
+                                                       pressure,
+                                                       temperature,
+                                                       compDens,
+                                                       pressureScalingFactor,
+                                                       compDensScalingFactor,
+                                                       temperatureScalingFactor,
+                                                       dofManager.rankOffset(),
+                                                       m_numComponents,
+                                                       dofKey,
+                                                       subRegion,
+                                                       localSolution,
+                                                       temperatureOffset )
+    : isothermalCompositionalMultiphaseBaseKernels::
+            SolutionScalingKernelFactory::
+            createAndLaunch< parallelDevicePolicy<> >( m_maxRelativePresChange,
+                                                       m_maxAbsolutePresChange,
+                                                       m_maxCompFracChange,
+                                                       m_maxRelativeCompDensChange,
+                                                       pressure,
+                                                       compDens,
+                                                       pressureScalingFactor,
+                                                       compDensScalingFactor,
+                                                       dofManager.rankOffset(),
+                                                       m_numComponents,
+                                                       dofKey,
+                                                       subRegion,
+                                                       localSolution );
+        if( subRegion.size() > 0 || subRegion.size() !=  subRegion.getNumberOfGhosts() )
+        {
+          if( m_scalingType == ScalingType::Global )
+          {
+            scalingFactor = std::min( scalingFactor, subRegionData.localMinVal );
+          }
+
+          regionDeltaPresMaxLoc.push_back( valueAndLocationType( subRegionData.localMaxDeltaPres, localToGlobalMap[subRegionData.localMaxDeltaPresLoc] ) );
+          minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.localMinPresScalingFactor );
+
+          regionDeltaCompDensMaxLoc.push_back( valueAndLocationType( subRegionData.localMaxDeltaCompDens, localToGlobalMap[subRegionData.localMaxDeltaCompDensLoc] ) );
+          minCompDensScalingFactor = std::min( minCompDensScalingFactor, subRegionData.localMinCompDensScalingFactor );
+
+          if( m_isThermal )
+          {
+            regionDeltaTempMaxLoc.push_back( valueAndLocationType( subRegionData.localMaxDeltaTemp, localToGlobalMap[subRegionData.localMaxDeltaTempLoc] ) );
+            minTempScalingFactor = std::min( minTempScalingFactor, subRegionData.localMinTempScalingFactor );
+          }
+        }
+      } );
+    } );
+
+    auto [localDeltaPresMax, localPresMaxLoc] = *std::max_element( begin( regionDeltaPresMaxLoc ), end( regionDeltaPresMaxLoc ), []( auto & lhs, auto & rhs )   {
+      return lhs.value   < rhs.value;
+    } );
+    auto globalDeltaPresMax = MpiWrapper::maxValLoc( valueAndLocationType( localDeltaPresMax, localPresMaxLoc ));
+    auto [ localDeltaCompDensMax, localCompDensMaxLoc ] = *std::max_element( begin( regionDeltaCompDensMaxLoc ), end( regionDeltaCompDensMaxLoc ), []( auto & lhs, auto & rhs )   {
+      return lhs.value   < rhs.value;
+    } );
+    auto globalDeltaCompDensMax = MpiWrapper::maxValLoc( valueAndLocationType( localDeltaCompDensMax, localCompDensMaxLoc ));
+
+    scalingFactor = MpiWrapper::min( scalingFactor );
+    minPresScalingFactor = MpiWrapper::min( minPresScalingFactor );
+    minCompDensScalingFactor = MpiWrapper::min( minCompDensScalingFactor );
 
+    string const massUnit = m_useMass ? "kg/m3" : "mol/m3";
+    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution,
+                                GEOS_FMT( "        {}: Max pressure change = {:.3f} Pa (before scaling) at cell {}",
+                                          getName(),
+                                          globalDeltaPresMax.value,
+                                          globalDeltaPresMax.location ) );
 
-  std::vector< valueAndLocationType > regionDeltaPresMaxLoc;
-  std::vector< valueAndLocationType > regionDeltaTempMaxLoc;
-  std::vector< valueAndLocationType > regionDeltaCompDensMaxLoc;
+    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution,
+                                GEOS_FMT( "        {}: Max component density change = {:.3f} {} (before scaling) at cell {}",
+                                          getName(),
+                                          globalDeltaCompDensMax.value,
+                                          massUnit,
+                                          globalDeltaCompDensMax.location ) );
+
+    if( m_isThermal )
+    {
+      auto [localDeltaTempMax, localDeltaTempMaxLoc  ] = *std::max_element( begin( regionDeltaTempMaxLoc ), end( regionDeltaTempMaxLoc ), []( auto & lhs, auto & rhs )   {
+        return lhs.value   < rhs.value;
+      } );
+      auto globalMaxDeltaTemp = MpiWrapper::maxValLoc( valueAndLocationType( localDeltaTempMax, localDeltaTempMaxLoc ));
+
+      minTempScalingFactor = MpiWrapper::min( minTempScalingFactor );
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution,
+                                  GEOS_FMT( "        {}: Max temperature change = {:.3f} K (before scaling) at cell maxRegionDeltaTempLoc {}",
+                                            getName(),
+                                            globalMaxDeltaTemp.value,
+                                            globalMaxDeltaTemp.location ) );
+    }
+
+    if( m_scalingType == ScalingType::Local )
+    {
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Min pressure scaling factor = {}", getName(), minPresScalingFactor ) );
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Min component density scaling factor = {}", getName(), minCompDensScalingFactor ) );
+      if( m_isThermal )
+      {
+        GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Min temperature scaling factor = {}", getName(), minTempScalingFactor ) );
+      }
+    }
+
+    return LvArray::math::max( scalingFactor, m_minScalingFactor );
+  }
+}
+
+real64 CompositionalMultiphaseFVM::scalingForSystemSolutionZFormulation( DomainPartition & domain,
+                                                                         DofManager const & dofManager,
+                                                                         arrayView1d< real64 const > const & localSolution )
+{
+  string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
+  real64 scalingFactor = 1.0;
+  real64 maxDeltaPres = 0.0, maxDeltaCompFrac = 0.0, maxDeltaTemp = 0.0;
+  real64 minPresScalingFactor = 1.0, minCompFracScalingFactor = 1.0, minTempScalingFactor = 1.0;
 
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
                                                                MeshLevel & mesh,
@@ -499,120 +685,65 @@ real64 CompositionalMultiphaseFVM::scalingForSystemSolution( DomainPartition & d
                                                 [&]( localIndex const,
                                                      ElementSubRegionBase & subRegion )
     {
-      arrayView1d< globalIndex const > const localToGlobalMap = subRegion.localToGlobalMap();
       arrayView1d< real64 const > const pressure = subRegion.getField< fields::flow::pressure >();
-      arrayView1d< real64 const > const temperature = subRegion.getField< fields::flow::temperature >();
-      arrayView2d< real64 const, compflow::USD_COMP > const compDens = subRegion.getField< fields::flow::globalCompDensity >();
+      //arrayView1d< real64 const > const temperature = subRegion.getField< fields::flow::temperature >();
+      arrayView2d< real64 const, compflow::USD_COMP > const compFrac = subRegion.getField< fields::flow::globalCompFraction >();
       arrayView1d< real64 > pressureScalingFactor = subRegion.getField< fields::flow::pressureScalingFactor >();
-      arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< fields::flow::temperatureScalingFactor >();
-      arrayView1d< real64 > compDensScalingFactor = subRegion.getField< fields::flow::globalCompDensityScalingFactor >();
-
-      const integer temperatureOffset = m_numComponents+1;
-
-      auto const subRegionData =
-        m_isThermal
-  ? thermalCompositionalMultiphaseBaseKernels::
-          SolutionScalingKernelFactory::
-          createAndLaunch< parallelDevicePolicy<> >( m_maxRelativePresChange,
-                                                     m_maxAbsolutePresChange,
-                                                     m_maxRelativeTempChange,
-                                                     m_maxCompFracChange,
-                                                     m_maxRelativeCompDensChange,
-                                                     pressure,
-                                                     temperature,
-                                                     compDens,
-                                                     pressureScalingFactor,
-                                                     compDensScalingFactor,
-                                                     temperatureScalingFactor,
-                                                     dofManager.rankOffset(),
-                                                     m_numComponents,
-                                                     dofKey,
-                                                     subRegion,
-                                                     localSolution,
-                                                     temperatureOffset )
-  : isothermalCompositionalMultiphaseBaseKernels::
-          SolutionScalingKernelFactory::
-          createAndLaunch< parallelDevicePolicy<> >( m_maxRelativePresChange,
-                                                     m_maxAbsolutePresChange,
-                                                     m_maxCompFracChange,
-                                                     m_maxRelativeCompDensChange,
-                                                     pressure,
-                                                     compDens,
-                                                     pressureScalingFactor,
-                                                     compDensScalingFactor,
-                                                     dofManager.rankOffset(),
-                                                     m_numComponents,
-                                                     dofKey,
-                                                     subRegion,
-                                                     localSolution );
-      if( subRegion.size() > 0 || subRegion.size() !=  subRegion.getNumberOfGhosts() )
+      //arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< fields::flow::temperatureScalingFactor >();
+      arrayView1d< real64 > compFracScalingFactor = subRegion.getField< fields::flow::globalCompFractionScalingFactor >();
+
+      auto const subRegionData = isothermalCompositionalMultiphaseBaseKernels::
+                                   SolutionScalingZFormulationKernelFactory::
+                                   createAndLaunch< parallelDevicePolicy<> >( m_maxRelativePresChange,
+                                                                              m_maxAbsolutePresChange,
+                                                                              m_maxCompFracChange,
+                                                                              pressure,
+                                                                              compFrac,
+                                                                              pressureScalingFactor,
+                                                                              compFracScalingFactor,
+                                                                              dofManager.rankOffset(),
+                                                                              m_numComponents,
+                                                                              dofKey,
+                                                                              subRegion,
+                                                                              localSolution );
+
+      if( m_scalingType == ScalingType::Global )
       {
-        if( m_scalingType == ScalingType::Global )
-        {
-          scalingFactor = std::min( scalingFactor, subRegionData.localMinVal );
-        }
-
-        regionDeltaPresMaxLoc.push_back( valueAndLocationType( subRegionData.localMaxDeltaPres, localToGlobalMap[subRegionData.localMaxDeltaPresLoc] ) );
-        minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.localMinPresScalingFactor );
-
-        regionDeltaCompDensMaxLoc.push_back( valueAndLocationType( subRegionData.localMaxDeltaCompDens, localToGlobalMap[subRegionData.localMaxDeltaCompDensLoc] ) );
-        minCompDensScalingFactor = std::min( minCompDensScalingFactor, subRegionData.localMinCompDensScalingFactor );
-
-        if( m_isThermal )
-        {
-          regionDeltaTempMaxLoc.push_back( valueAndLocationType( subRegionData.localMaxDeltaTemp, localToGlobalMap[subRegionData.localMaxDeltaTempLoc] ) );
-          minTempScalingFactor = std::min( minTempScalingFactor, subRegionData.localMinTempScalingFactor );
-        }
+        scalingFactor = std::min( scalingFactor, subRegionData.localMinVal );
       }
-    } );
-  } );
+      maxDeltaPres  = std::max( maxDeltaPres, subRegionData.localMaxDeltaPres );
+      maxDeltaCompFrac = std::max( maxDeltaCompFrac, subRegionData.localMaxDeltaCompFrac );
+      maxDeltaTemp = std::max( maxDeltaTemp, subRegionData.localMaxDeltaTemp );
+      minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.localMinPresScalingFactor );
+      minCompFracScalingFactor = std::min( minCompFracScalingFactor, subRegionData.localMinCompFracScalingFactor );
+      minTempScalingFactor = std::min( minTempScalingFactor, subRegionData.localMinTempScalingFactor );
 
-  auto [localDeltaPresMax, localPresMaxLoc] = *std::max_element( begin( regionDeltaPresMaxLoc ), end( regionDeltaPresMaxLoc ), []( auto & lhs, auto & rhs )   {
-    return lhs.value   < rhs.value;
-  } );
-  auto globalDeltaPresMax = MpiWrapper::maxValLoc( valueAndLocationType( localDeltaPresMax, localPresMaxLoc ));
-  auto [ localDeltaCompDensMax, localCompDensMaxLoc ] = *std::max_element( begin( regionDeltaCompDensMaxLoc ), end( regionDeltaCompDensMaxLoc ), []( auto & lhs, auto & rhs )   {
-    return lhs.value   < rhs.value;
+    } );
   } );
-  auto globalDeltaCompDensMax = MpiWrapper::maxValLoc( valueAndLocationType( localDeltaCompDensMax, localCompDensMaxLoc ));
 
   scalingFactor = MpiWrapper::min( scalingFactor );
+  maxDeltaPres  = MpiWrapper::max( maxDeltaPres );
+  maxDeltaCompFrac = MpiWrapper::max( maxDeltaCompFrac );
   minPresScalingFactor = MpiWrapper::min( minPresScalingFactor );
-  minCompDensScalingFactor = MpiWrapper::min( minCompDensScalingFactor );
-
-  string const massUnit = m_useMass ? "kg/m3" : "mol/m3";
-  GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution,
-                              GEOS_FMT( "        {}: Max pressure change = {:.3f} Pa (before scaling) at cell {}",
-                                        getName(),
-                                        globalDeltaPresMax.value,
-                                        globalDeltaPresMax.location ) );
-
-  GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution,
-                              GEOS_FMT( "        {}: Max component density change = {:.3f} {} (before scaling) at cell {}",
-                                        getName(),
-                                        globalDeltaCompDensMax.value,
-                                        massUnit,
-                                        globalDeltaCompDensMax.location ) );
+  minCompFracScalingFactor = MpiWrapper::min( minCompFracScalingFactor );
+
+  GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Max pressure change = {} Pa (before scaling)",
+                                                           getName(), GEOS_FMT( "{:.{}f}", maxDeltaPres, 3 ) ) );
+  GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Max component fraction change = {} (before scaling)",
+                                                           getName(), GEOS_FMT( "{:.{}f}", maxDeltaCompFrac, 3 ) ) );
 
   if( m_isThermal )
   {
-    auto [localDeltaTempMax, localDeltaTempMaxLoc  ] = *std::max_element( begin( regionDeltaTempMaxLoc ), end( regionDeltaTempMaxLoc ), []( auto & lhs, auto & rhs )   {
-      return lhs.value   < rhs.value;
-    } );
-    auto globalMaxDeltaTemp = MpiWrapper::maxValLoc( valueAndLocationType( localDeltaTempMax, localDeltaTempMaxLoc ));
-
+    maxDeltaTemp = MpiWrapper::max( maxDeltaTemp );
     minTempScalingFactor = MpiWrapper::min( minTempScalingFactor );
-    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution,
-                                GEOS_FMT( "        {}: Max temperature change = {:.3f} K (before scaling) at cell maxRegionDeltaTempLoc {}",
-                                          getName(),
-                                          globalMaxDeltaTemp.value,
-                                          globalMaxDeltaTemp.location ) );
+    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Max temperature change = {} K (before scaling)",
+                                                             getName(), GEOS_FMT( "{:.{}f}", maxDeltaTemp, 3 ) ) );
   }
 
   if( m_scalingType == ScalingType::Local )
   {
     GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Min pressure scaling factor = {}", getName(), minPresScalingFactor ) );
-    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Min component density scaling factor = {}", getName(), minCompDensScalingFactor ) );
+    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Min component fraction scaling factor = {}", getName(), minCompFracScalingFactor ) );
     if( m_isThermal )
     {
       GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Min temperature scaling factor = {}", getName(), minTempScalingFactor ) );
@@ -629,97 +760,105 @@ bool CompositionalMultiphaseFVM::checkSystemSolution( DomainPartition & domain,
 {
   GEOS_MARK_FUNCTION;
 
-  string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
-  integer localCheck = 1;
-  real64 minPres = 0.0, minDens = 0.0, minTotalDens = 0.0;
-  integer numNegPres = 0, numNegDens = 0, numNegTotalDens = 0;
-
-  forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
-                                                               MeshLevel & mesh,
-                                                               arrayView1d< string const > const & regionNames )
+  // TO DO: Implement the solution check for Z Formulation
+  if( m_useZFormulation )
   {
-    mesh.getElemManager().forElementSubRegions( regionNames,
-                                                [&]( localIndex const,
-                                                     ElementSubRegionBase & subRegion )
+    return true;
+  }
+  else
+  {
+    string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
+    integer localCheck = 1;
+    real64 minPres = 0.0, minDens = 0.0, minTotalDens = 0.0;
+    integer numNegPres = 0, numNegDens = 0, numNegTotalDens = 0;
+
+    forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+                                                                 MeshLevel & mesh,
+                                                                 arrayView1d< string const > const & regionNames )
     {
-      arrayView1d< real64 const > const pressure =
-        subRegion.getField< fields::flow::pressure >();
-      arrayView1d< real64 const > const temperature =
-        subRegion.getField< fields::flow::temperature >();
-      arrayView2d< real64 const, compflow::USD_COMP > const compDens =
-        subRegion.getField< fields::flow::globalCompDensity >();
-      arrayView1d< real64 > pressureScalingFactor = subRegion.getField< fields::flow::pressureScalingFactor >();
-      arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< fields::flow::temperatureScalingFactor >();
-      arrayView1d< real64 > compDensScalingFactor = subRegion.getField< fields::flow::globalCompDensityScalingFactor >();
-      // check that pressure and component densities are non-negative
-      // for thermal, check that temperature is above 273.15 K
-      const integer temperatureOffset = m_numComponents+1;
-      auto const subRegionData =
-        m_isThermal
-  ? thermalCompositionalMultiphaseBaseKernels::
-          SolutionCheckKernelFactory::
-          createAndLaunch< parallelDevicePolicy<> >( m_allowCompDensChopping,
-                                                     m_allowNegativePressure,
-                                                     m_scalingType,
-                                                     scalingFactor,
-                                                     pressure,
-                                                     temperature,
-                                                     compDens,
-                                                     pressureScalingFactor,
-                                                     temperatureScalingFactor,
-                                                     compDensScalingFactor,
-                                                     dofManager.rankOffset(),
-                                                     m_numComponents,
-                                                     dofKey,
-                                                     subRegion,
-                                                     localSolution,
-                                                     temperatureOffset )
-  : isothermalCompositionalMultiphaseBaseKernels::
-          SolutionCheckKernelFactory::
-          createAndLaunch< parallelDevicePolicy<> >( m_allowCompDensChopping,
-                                                     m_allowNegativePressure,
-                                                     m_scalingType,
-                                                     scalingFactor,
-                                                     pressure,
-                                                     compDens,
-                                                     pressureScalingFactor,
-                                                     compDensScalingFactor,
-                                                     dofManager.rankOffset(),
-                                                     m_numComponents,
-                                                     dofKey,
-                                                     subRegion,
-                                                     localSolution );
-
-      localCheck = std::min( localCheck, subRegionData.localMinVal );
-
-      minPres  = std::min( minPres, subRegionData.localMinPres );
-      minDens = std::min( minDens, subRegionData.localMinDens );
-      minTotalDens = std::min( minTotalDens, subRegionData.localMinTotalDens );
-      numNegPres += subRegionData.localNumNegPressures;
-      numNegDens += subRegionData.localNumNegDens;
-      numNegTotalDens += subRegionData.localNumNegTotalDens;
+      mesh.getElemManager().forElementSubRegions( regionNames,
+                                                  [&]( localIndex const,
+                                                       ElementSubRegionBase & subRegion )
+      {
+        arrayView1d< real64 const > const pressure =
+          subRegion.getField< fields::flow::pressure >();
+        arrayView1d< real64 const > const temperature =
+          subRegion.getField< fields::flow::temperature >();
+        arrayView2d< real64 const, compflow::USD_COMP > const compDens =
+          subRegion.getField< fields::flow::globalCompDensity >();
+        arrayView1d< real64 > pressureScalingFactor = subRegion.getField< fields::flow::pressureScalingFactor >();
+        arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< fields::flow::temperatureScalingFactor >();
+        arrayView1d< real64 > compDensScalingFactor = subRegion.getField< fields::flow::globalCompDensityScalingFactor >();
+        // check that pressure and component densities are non-negative
+        // for thermal, check that temperature is above 273.15 K
+        const integer temperatureOffset = m_numComponents+1;
+        auto const subRegionData =
+          m_isThermal
+    ? thermalCompositionalMultiphaseBaseKernels::
+            SolutionCheckKernelFactory::
+            createAndLaunch< parallelDevicePolicy<> >( m_allowCompDensChopping,
+                                                       m_allowNegativePressure,
+                                                       m_scalingType,
+                                                       scalingFactor,
+                                                       pressure,
+                                                       temperature,
+                                                       compDens,
+                                                       pressureScalingFactor,
+                                                       temperatureScalingFactor,
+                                                       compDensScalingFactor,
+                                                       dofManager.rankOffset(),
+                                                       m_numComponents,
+                                                       dofKey,
+                                                       subRegion,
+                                                       localSolution,
+                                                       temperatureOffset )
+    : isothermalCompositionalMultiphaseBaseKernels::
+            SolutionCheckKernelFactory::
+            createAndLaunch< parallelDevicePolicy<> >( m_allowCompDensChopping,
+                                                       m_allowNegativePressure,
+                                                       m_scalingType,
+                                                       scalingFactor,
+                                                       pressure,
+                                                       compDens,
+                                                       pressureScalingFactor,
+                                                       compDensScalingFactor,
+                                                       dofManager.rankOffset(),
+                                                       m_numComponents,
+                                                       dofKey,
+                                                       subRegion,
+                                                       localSolution );
+
+        localCheck = std::min( localCheck, subRegionData.localMinVal );
+
+        minPres  = std::min( minPres, subRegionData.localMinPres );
+        minDens = std::min( minDens, subRegionData.localMinDens );
+        minTotalDens = std::min( minTotalDens, subRegionData.localMinTotalDens );
+        numNegPres += subRegionData.localNumNegPressures;
+        numNegDens += subRegionData.localNumNegDens;
+        numNegTotalDens += subRegionData.localNumNegTotalDens;
+      } );
     } );
-  } );
 
-  minPres  = MpiWrapper::min( minPres );
-  minDens = MpiWrapper::min( minDens );
-  minTotalDens = MpiWrapper::min( minTotalDens );
-  numNegPres = MpiWrapper::sum( numNegPres );
-  numNegDens = MpiWrapper::sum( numNegDens );
-  numNegTotalDens = MpiWrapper::sum( numNegTotalDens );
-
-  if( numNegPres > 0 )
-    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Number of negative pressure values: {}, minimum value: {} Pa",
-                                                             getName(), numNegPres, fmt::format( "{:.{}f}", minPres, 3 ) ) );
-  string const massUnit = m_useMass ? "kg/m3" : "mol/m3";
-  if( numNegDens > 0 )
-    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Number of negative component density values: {}, minimum value: {} {}}",
-                                                             getName(), numNegDens, fmt::format( "{:.{}f}", minDens, 3 ), massUnit ) );
-  if( minTotalDens > 0 )
-    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Number of negative total density values: {}, minimum value: {} {}}",
-                                                             getName(), minTotalDens, fmt::format( "{:.{}f}", minDens, 3 ), massUnit ) );
-
-  return MpiWrapper::min( localCheck );
+    minPres  = MpiWrapper::min( minPres );
+    minDens = MpiWrapper::min( minDens );
+    minTotalDens = MpiWrapper::min( minTotalDens );
+    numNegPres = MpiWrapper::sum( numNegPres );
+    numNegDens = MpiWrapper::sum( numNegDens );
+    numNegTotalDens = MpiWrapper::sum( numNegTotalDens );
+
+    if( numNegPres > 0 )
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Number of negative pressure values: {}, minimum value: {} Pa",
+                                                               getName(), numNegPres, fmt::format( "{:.{}f}", minPres, 3 ) ) );
+    string const massUnit = m_useMass ? "kg/m3" : "mol/m3";
+    if( numNegDens > 0 )
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Number of negative component density values: {}, minimum value: {} {}}",
+                                                               getName(), numNegDens, fmt::format( "{:.{}f}", minDens, 3 ), massUnit ) );
+    if( minTotalDens > 0 )
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Number of negative total density values: {}, minimum value: {} {}}",
+                                                               getName(), minTotalDens, fmt::format( "{:.{}f}", minDens, 3 ), massUnit ) );
+
+    return MpiWrapper::min( localCheck );
+  }
 }
 
 void CompositionalMultiphaseFVM::applySystemSolution( DofManager const & dofManager,
@@ -752,21 +891,43 @@ void CompositionalMultiphaseFVM::applySystemSolution( DofManager const & dofMana
                                  pressureMask );
   }
 
-  if( localScaling )
+  if( m_useZFormulation )
   {
-    dofManager.addVectorToField( localSolution,
-                                 viewKeyStruct::elemDofFieldString(),
-                                 fields::flow::globalCompDensity::key(),
-                                 fields::flow::globalCompDensityScalingFactor::key(),
-                                 componentMask );
+    if( localScaling )
+    {
+      dofManager.addVectorToField( localSolution,
+                                   viewKeyStruct::elemDofFieldString(),
+                                   fields::flow::globalCompFraction::key(),
+                                   fields::flow::globalCompFractionScalingFactor::key(),
+                                   componentMask );
+    }
+    else
+    {
+      dofManager.addVectorToField( localSolution,
+                                   viewKeyStruct::elemDofFieldString(),
+                                   fields::flow::globalCompFraction::key(),
+                                   scalingFactor,
+                                   componentMask );
+    }
   }
   else
   {
-    dofManager.addVectorToField( localSolution,
-                                 viewKeyStruct::elemDofFieldString(),
-                                 fields::flow::globalCompDensity::key(),
-                                 scalingFactor,
-                                 componentMask );
+    if( localScaling )
+    {
+      dofManager.addVectorToField( localSolution,
+                                   viewKeyStruct::elemDofFieldString(),
+                                   fields::flow::globalCompDensity::key(),
+                                   fields::flow::globalCompDensityScalingFactor::key(),
+                                   componentMask );
+    }
+    else
+    {
+      dofManager.addVectorToField( localSolution,
+                                   viewKeyStruct::elemDofFieldString(),
+                                   fields::flow::globalCompDensity::key(),
+                                   scalingFactor,
+                                   componentMask );
+    }
   }
 
   if( m_isThermal )
@@ -792,16 +953,21 @@ 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
+
   if( m_allowCompDensChopping )
   {
-    chopNegativeDensities( domain );
+    if( m_useZFormulation )
+      chopNegativeCompFractions( domain );
+    else
+      chopNegativeDensities( domain );
   }
 
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
                                                                MeshLevel & mesh,
                                                                arrayView1d< string const > const & regionNames )
   {
-    std::vector< string > fields{ fields::flow::pressure::key(), 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() );
@@ -824,10 +990,11 @@ void CompositionalMultiphaseFVM::updatePhaseMobility( ObjectManagerBase & dataGr
   string const & relpermName = dataGroup.getReference< string >( viewKeyStruct::relPermNamesString() );
   RelativePermeabilityBase const & relperm = getConstitutiveModel< RelativePermeabilityBase >( dataGroup, relpermName );
 
-  if( m_isThermal )
+  if( m_useZFormulation )
   {
-    thermalCompositionalMultiphaseFVMKernels::
-      PhaseMobilityKernelFactory::
+    // For now: isothermal only
+    isothermalCompositionalMultiphaseFVMKernels::
+      PhaseMobilityZFormulationKernelFactory::
       createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
                                                  m_numPhases,
                                                  dataGroup,
@@ -836,13 +1003,26 @@ void CompositionalMultiphaseFVM::updatePhaseMobility( ObjectManagerBase & dataGr
   }
   else
   {
-    isothermalCompositionalMultiphaseFVMKernels::
-      PhaseMobilityKernelFactory::
-      createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
-                                                 m_numPhases,
-                                                 dataGroup,
-                                                 fluid,
-                                                 relperm );
+    if( m_isThermal )
+    {
+      thermalCompositionalMultiphaseFVMKernels::
+        PhaseMobilityKernelFactory::
+        createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                   m_numPhases,
+                                                   dataGroup,
+                                                   fluid,
+                                                   relperm );
+    }
+    else
+    {
+      isothermalCompositionalMultiphaseFVMKernels::
+        PhaseMobilityKernelFactory::
+        createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                   m_numPhases,
+                                                   dataGroup,
+                                                   fluid,
+                                                   relperm );
+    }
   }
 }
 
@@ -1005,7 +1185,7 @@ void CompositionalMultiphaseFVM::applyFaceDirichletBC( real64 const time_n,
   if( m_nonlinearSolverParameters.m_numNewtonIterations == 0 )
   {
     bool const bcConsistent = validateFaceDirichletBC( domain, time_n + dt );
-    GEOS_ERROR_IF( !bcConsistent, GEOS_FMT( "CompositionalMultiphaseBase {}: inconsistent boundary conditions", getDataContext() ) );
+    GEOS_ERROR_IF( !bcConsistent, GEOS_FMT( "{}: inconsistent boundary conditions", getDataContext() ) );
   }
 
   FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance();
@@ -1062,11 +1242,10 @@ void CompositionalMultiphaseFVM::applyFaceDirichletBC( real64 const time_n,
 
       string const & elemDofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
 
-      if( m_isThermal )
+      if( m_useZFormulation )
       {
-        //todo (jafranc) extend upwindScheme name if satisfied in isothermalCase
-        thermalCompositionalMultiphaseFVMKernels::
-          DirichletFluxComputeKernelFactory::
+        isothermalCompositionalMultiphaseFVMKernels::
+          DirichletFluxComputeZFormulationKernelFactory::
           createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
                                                      m_numPhases,
                                                      dofManager.rankOffset(),
@@ -1080,24 +1259,47 @@ void CompositionalMultiphaseFVM::applyFaceDirichletBC( real64 const time_n,
                                                      dt,
                                                      localMatrix,
                                                      localRhs );
+
       }
       else
       {
-        isothermalCompositionalMultiphaseFVMKernels::
-          DirichletFluxComputeKernelFactory::
-          createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
-                                                     m_numPhases,
-                                                     dofManager.rankOffset(),
-                                                     m_useTotalMassEquation,
-                                                     elemDofKey,
-                                                     getName(),
-                                                     faceManager,
-                                                     elemManager,
-                                                     stencilWrapper,
-                                                     multiFluidBase,
-                                                     dt,
-                                                     localMatrix,
-                                                     localRhs );
+        if( m_isThermal )
+        {
+          //todo (jafranc) extend upwindScheme name if satisfied in isothermalCase
+          thermalCompositionalMultiphaseFVMKernels::
+            DirichletFluxComputeKernelFactory::
+            createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                       m_numPhases,
+                                                       dofManager.rankOffset(),
+                                                       m_useTotalMassEquation,
+                                                       elemDofKey,
+                                                       getName(),
+                                                       faceManager,
+                                                       elemManager,
+                                                       stencilWrapper,
+                                                       multiFluidBase,
+                                                       dt,
+                                                       localMatrix,
+                                                       localRhs );
+        }
+        else
+        {
+          isothermalCompositionalMultiphaseFVMKernels::
+            DirichletFluxComputeKernelFactory::
+            createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
+                                                       m_numPhases,
+                                                       dofManager.rankOffset(),
+                                                       m_useTotalMassEquation,
+                                                       elemDofKey,
+                                                       getName(),
+                                                       faceManager,
+                                                       elemManager,
+                                                       stencilWrapper,
+                                                       multiFluidBase,
+                                                       dt,
+                                                       localMatrix,
+                                                       localRhs );
+        }
       }
 
     } );
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp
index f14b243674..f1437bd45f 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp
@@ -109,6 +109,11 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase
                             DofManager const & dofManager,
                             arrayView1d< real64 const > const & localSolution ) override;
 
+  real64
+  scalingForSystemSolutionZFormulation( DomainPartition & domain,
+                                        DofManager const & dofManager,
+                                        arrayView1d< real64 const > const & localSolution );
+
   virtual bool
   checkSystemSolution( DomainPartition & domain,
                        DofManager const & dofManager,
@@ -138,7 +143,6 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase
                                CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                arrayView1d< real64 > const & localRhs ) const override;
 
-
   virtual void
   updatePhaseMobility( ObjectManagerBase & dataGroup ) const override;
 
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp
index 021ce8fd96..922b8fbf01 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp
@@ -509,7 +509,6 @@ real64 CompositionalMultiphaseHybridFVM::scalingForSystemSolution( DomainPartiti
   return LvArray::math::max( MpiWrapper::min( scalingFactor ), m_minScalingFactor );
 }
 
-
 bool CompositionalMultiphaseHybridFVM::checkSystemSolution( DomainPartition & domain,
                                                             DofManager const & dofManager,
                                                             arrayView1d< real64 const > const & localSolution,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBLFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBLFields.hpp
index d83c199f8b..7af784da50 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBLFields.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBLFields.hpp
@@ -38,14 +38,6 @@ using array2dLayoutComp = array2d< real64, compflow::LAYOUT_COMP >;
 using array2dLayoutOBLOpVals = array2d< real64, compflow::LAYOUT_OBL_OPERATOR_VALUES >;
 using array3dLayoutOBLOpDers = array3d< real64, compflow::LAYOUT_OBL_OPERATOR_DERIVATIVES >;
 
-DECLARE_FIELD( globalCompFraction_n,
-               "globalCompFraction_n",
-               array2dLayoutComp,
-               0,
-               NOPLOT,
-               NO_WRITE,
-               "Global component fraction at the previous converged time step" );
-
 DECLARE_FIELD( bcGlobalCompFraction,
                "bcGlobalCompFraction",
                array2dLayoutComp,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/AccumulationZFormulationKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/AccumulationZFormulationKernel.hpp
new file mode 100644
index 0000000000..29029f6a7d
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/AccumulationZFormulationKernel.hpp
@@ -0,0 +1,417 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file AccumulationZFormulationKernel.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_ACCUMULATIONZFORMULATIONKERNEL_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_ACCUMULATIONZFORMULATIONKERNEL_HPP
+
+#include "codingUtilities/Utilities.hpp"
+#include "common/DataLayouts.hpp"
+#include "common/DataTypes.hpp"
+#include "common/GEOS_RAJA_Interface.hpp"
+#include "constitutive/solid/CoupledSolidBase.hpp"
+#include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
+#include "mesh/ElementSubRegionBase.hpp"
+#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
+#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp"
+#include "physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/KernelLaunchSelectors.hpp"
+
+namespace geos
+{
+
+namespace isothermalCompositionalMultiphaseBaseKernels
+{
+
+static constexpr real64 minCompFracForDivision = 0;
+
+/******************************** AccumulationKernel ********************************/
+
+/**
+ * @class AccumulationZFormulationKernel
+ * @tparam NUM_COMP number of fluid components
+ * @tparam NUM_DOF number of degrees of freedom
+ * @brief Define the interface for the assembly kernel in charge of accumulation and volume balance
+ */
+template< integer NUM_COMP, integer NUM_DOF >
+class AccumulationZFormulationKernel
+{
+public:
+
+  /// Compile time value for the number of components
+  static constexpr integer numComp = NUM_COMP;
+
+  /// Compute time value for the number of degrees of freedom
+  static constexpr integer numDof = NUM_DOF;
+
+  /// Compute time value for the number of equations
+  static constexpr integer numEqn = NUM_DOF;
+
+  /**
+   * @brief Constructor
+   * @param[in] numPhases the number of fluid phases
+   * @param[in] rankOffset the offset of my MPI rank
+   * @param[in] dofKey the string key to retrieve the degress of freedom numbers
+   * @param[in] subRegion the element subregion
+   * @param[in] fluid the fluid model
+   * @param[in] solid the solid model
+   * @param[inout] localMatrix the local CRS matrix
+   * @param[inout] localRhs the local right-hand side vector
+   */
+  AccumulationZFormulationKernel( localIndex const numPhases,
+                                  globalIndex const rankOffset,
+                                  string const dofKey,
+                                  ElementSubRegionBase const & subRegion,
+                                  constitutive::MultiFluidBase const & fluid,
+                                  constitutive::CoupledSolidBase const & solid,
+                                  CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                                  arrayView1d< real64 > const & localRhs,
+                                  BitFlags< KernelFlags > const KernelFlags )
+    : m_numPhases( numPhases ),
+    m_rankOffset( rankOffset ),
+    m_dofNumber( subRegion.getReference< array1d< globalIndex > >( dofKey ) ),
+    m_elemGhostRank( subRegion.ghostRank() ),
+    m_volume( subRegion.getElementVolume() ),
+    m_porosity( solid.getPorosity() ),
+    m_dPoro_dPres( solid.getDporosity_dPressure() ),
+    m_phaseVolFrac( subRegion.getField< fields::flow::phaseVolumeFraction >() ),
+    m_dPhaseVolFrac( subRegion.getField< fields::flow::dPhaseVolumeFraction >() ),
+    m_totalDens( fluid.totalDensity() ),
+    m_dTotalDens( fluid.dTotalDensity() ),
+    m_compFrac( subRegion.getField< fields::flow::globalCompFraction >() ),
+    m_compAmount_n( subRegion.getField< fields::flow::compAmount_n >() ),
+    m_localMatrix( localMatrix ),
+    m_localRhs( localRhs ),
+    m_KernelFlags( KernelFlags )
+  {}
+
+  /**
+   * @struct StackVariables
+   * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack
+   */
+  struct StackVariables
+  {
+public:
+
+    // Pore volume information (used by both accumulation and volume balance)
+
+    /// Pore volume at time n+1
+    real64 poreVolume = 0.0;
+
+    /// Derivative of pore volume with respect to pressure
+    real64 dPoreVolume_dPres = 0.0;
+
+    // Residual information
+
+    /// Index of the local row corresponding to this element
+    localIndex localRow = -1;
+
+    /// Indices of the matrix rows/columns corresponding to the dofs in this element
+    globalIndex dofIndices[numDof]{};
+
+    /// C-array storage for the element local residual vector (all equations except volume balance)
+    real64 localResidual[numEqn]{};
+
+    /// C-array storage for the element local Jacobian matrix (all equations except volume balance, all dofs)
+    real64 localJacobian[numEqn][numDof]{};
+
+  };
+
+  /**
+   * @brief Getter for the ghost rank of an element
+   * @param[in] ei the element index
+   * @return the ghost rank of the element
+   */
+  GEOS_HOST_DEVICE
+  integer elemGhostRank( localIndex const ei ) const
+  { return m_elemGhostRank( ei ); }
+
+
+  /**
+   * @brief Performs the setup phase for the kernel.
+   * @param[in] ei the element index
+   * @param[in] stack the stack variables
+   */
+  GEOS_HOST_DEVICE
+  void setup( localIndex const ei,
+              StackVariables & stack ) const
+  {
+    // initialize the pore volume
+    stack.poreVolume = m_volume[ei] * m_porosity[ei][0];
+    stack.dPoreVolume_dPres = m_volume[ei] * m_dPoro_dPres[ei][0];
+
+    // set row index and degrees of freedom indices for this element
+    stack.localRow = m_dofNumber[ei] - m_rankOffset;
+    for( integer idof = 0; idof < numDof; ++idof )
+    {
+      stack.dofIndices[idof] = m_dofNumber[ei] + idof;
+    }
+  }
+
+  /**
+   * @brief Compute the local accumulation contributions to the residual and Jacobian
+   * @tparam FUNC the type of the function that can be used to customize the kernel
+   * @param[in] ei the element index
+   * @param[inout] stack the stack variables
+   * @param[in] phaseAmountKernelOp the function used to customize the kernel
+   */
+  template< typename FUNC = NoOpFunc >
+  GEOS_HOST_DEVICE
+  void computeAccumulation( localIndex const ei,
+                            StackVariables & stack ) const
+  {
+    using Deriv = constitutive::multifluid::DerivativeOffset;
+
+    arraySlice1d< real64 const, compflow::USD_COMP - 1 > const compFrac = m_compFrac[ei];
+    real64 const totalDensity = m_totalDens[ei][0];
+    arraySlice1d< real64 const, constitutive::multifluid::USD_FLUID_DC - 2 > const dTotalDens = m_dTotalDens[ei][0];
+
+    // ic - index of component whose conservation equation is assembled
+    // (i.e. row number in local matrix)
+    for( integer ic = 0; ic < numComp; ++ic )
+    {
+      real64 const compAmount = stack.poreVolume * totalDensity * compFrac[ic];
+      real64 const compAmount_n = m_compAmount_n[ei][ic];
+
+      stack.localResidual[ic] += compAmount - compAmount_n;
+
+      // derivatives with respect to pressure (p)
+      real64 const dCompAmount_dP = compFrac[ic] * (stack.dPoreVolume_dPres * totalDensity + stack.poreVolume * dTotalDens[Deriv::dP]);
+      stack.localJacobian[ic][0] += dCompAmount_dP;
+
+      // derivatives with respect to global component fraction (zc)
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        real64 dCompAmount_dC;
+        if( ic == jc )
+          dCompAmount_dC = stack.poreVolume * (totalDensity + dTotalDens[Deriv::dC+jc] * compFrac[ic]);
+        else
+          dCompAmount_dC = stack.poreVolume * (dTotalDens[Deriv::dC+jc] * compFrac[ic]);
+
+        stack.localJacobian[ic][jc + 1] += dCompAmount_dC;
+      }
+    }
+  }
+
+  /**
+   * @brief Compute the local volume balance contributions to the residual and Jacobian
+   * @tparam FUNC the type of the function that can be used to customize the kernel
+   * @param[in] ei the element index
+   * @param[inout] stack the stack variables
+   * @param[in] phaseVolFractionSumKernelOp the function used to customize the kernel
+   */
+  template< typename FUNC = NoOpFunc >
+  GEOS_HOST_DEVICE
+  void computeCompFracSum( localIndex const ei,
+                           StackVariables & stack ) const
+  {
+
+    arraySlice1d< real64 const, compflow::USD_PHASE - 1 > compFrac = m_compFrac[ei];
+
+    real64 oneMinusCompFracSum = 1.0;
+
+    // sum contributions to component accumulation from each component
+    for( integer ic = 0; ic < numComp; ++ic )
+    {
+      oneMinusCompFracSum -= compFrac[ic];
+      // no derivatives w.r.t pressure
+
+      // derivative w.r.t component fractions are unity: dzc_dzc = 1
+      stack.localJacobian[numComp][ic+1] -= 1;
+    }
+
+    //phaseVolFractionSumKernelOp( oneMinusCompFracSum );
+
+    // scale componentFraction-based volume balance by pore volume (for better scaling w.r.t. other equations)
+    stack.localResidual[numComp] = stack.poreVolume * oneMinusCompFracSum;
+    for( integer idof = 0; idof < numDof; ++idof )
+    {
+      stack.localJacobian[numComp][idof] *= stack.poreVolume;
+    }
+    stack.localJacobian[numComp][0] += stack.dPoreVolume_dPres * oneMinusCompFracSum;
+  }
+
+  /**
+   * @brief Performs the complete phase for the kernel.
+   * @param[in] ei the element index
+   * @param[inout] stack the stack variables
+   */
+  GEOS_HOST_DEVICE
+  void complete( localIndex const GEOS_UNUSED_PARAM( ei ),
+                 StackVariables & stack ) const
+  {
+    using namespace compositionalMultiphaseUtilities;
+
+    if( m_KernelFlags.isSet( KernelFlags::TotalMassEquation ) )
+    {
+      // apply equation/variable change transformation to the component mass balance equations
+      real64 work[numDof]{};
+      shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numDof, stack.localJacobian, work );
+      shiftElementsAheadByOneAndReplaceFirstElementWithSum( numComp, stack.localResidual );
+    }
+
+    // add contribution to residual and jacobian into:
+    // - the component mass balance equations (i = 0 to i = numComp-1)
+    // - the volume balance equations (i = numComp)
+    // note that numDof includes derivatives wrt temperature if this class is derived in ThermalKernels
+    integer const numRows = numComp+1;
+    for( integer i = 0; i < numRows; ++i )
+    {
+      m_localRhs[stack.localRow + i] += stack.localResidual[i];
+      m_localMatrix.addToRow< serialAtomic >( stack.localRow + i,
+                                              stack.dofIndices,
+                                              stack.localJacobian[i],
+                                              numDof );
+    }
+  }
+
+  /**
+   * @brief Performs the kernel launch
+   * @tparam POLICY the policy used in the RAJA kernels
+   * @tparam KERNEL_TYPE the kernel type
+   * @param[in] numElems the number of elements
+   * @param[inout] kernelComponent the kernel component providing access to setup/compute/complete functions and stack variables
+   */
+  template< typename POLICY, typename KERNEL_TYPE >
+  static void
+  launch( localIndex const numElems,
+          KERNEL_TYPE const & kernelComponent )
+  {
+    GEOS_MARK_FUNCTION;
+
+    forAll< POLICY >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei )
+    {
+      if( kernelComponent.elemGhostRank( ei ) >= 0 )
+      {
+        return;
+      }
+
+      typename KERNEL_TYPE::StackVariables stack;
+
+      kernelComponent.setup( ei, stack );
+      kernelComponent.computeAccumulation( ei, stack );
+      kernelComponent.computeCompFracSum( ei, stack );
+      kernelComponent.complete( ei, stack );
+    } );
+  }
+
+protected:
+
+  /// Number of fluid phases
+  integer const m_numPhases;
+
+  /// Offset for my MPI rank
+  globalIndex const m_rankOffset;
+
+  /// View on the dof numbers
+  arrayView1d< globalIndex const > const m_dofNumber;
+
+  /// View on the ghost ranks
+  arrayView1d< integer const > const m_elemGhostRank;
+
+  /// View on the element volumes
+  arrayView1d< real64 const > const m_volume;
+
+  /// Views on the porosity
+  arrayView2d< real64 const > const m_porosity;
+  arrayView2d< real64 const > const m_dPoro_dPres;
+
+  /// Views on the derivatives of comp fractions wrt component density
+  arrayView3d< real64 const, compflow::USD_COMP_DC > const m_dCompFrac_dCompDens;
+
+  /// Views on the phase volume fractions
+  arrayView2d< real64 const, compflow::USD_PHASE > const m_phaseVolFrac;
+  arrayView3d< real64 const, compflow::USD_PHASE_DC > const m_dPhaseVolFrac;
+
+  /// Views on the total density
+  arrayView2d< real64 const, constitutive::multifluid::USD_FLUID > const m_totalDens;
+  arrayView3d< real64 const, constitutive::multifluid::USD_FLUID_DC > const m_dTotalDens;
+
+  // View on component densities and component fractions
+  arrayView2d< real64 const, compflow::USD_COMP > m_compFrac;
+
+  // View on component amount (mass/moles) from previous time step
+  arrayView2d< real64 const, compflow::USD_COMP > m_compAmount_n;
+
+  /// View on the local CRS matrix
+  CRSMatrixView< real64, globalIndex const > const m_localMatrix;
+  /// View on the local RHS
+  arrayView1d< real64 > const m_localRhs;
+
+  BitFlags< KernelFlags > const m_KernelFlags;
+};
+
+/**
+ * @class AccumulationZFormulationKernelFactory
+ */
+class AccumulationZFormulationKernelFactory
+{
+public:
+
+  /**
+   * @brief Create a new kernel and launch
+   * @tparam POLICY the policy used in the RAJA kernel
+   * @param[in] numComps the number of fluid components
+   * @param[in] numPhases the number of fluid phases
+   * @param[in] rankOffset the offset of my MPI rank
+   * @param[in] dofKey the string key to retrieve the degress of freedom numbers
+   * @param[in] subRegion the element subregion
+   * @param[in] fluid the fluid model
+   * @param[in] solid the solid model
+   * @param[inout] localMatrix the local CRS matrix
+   * @param[inout] localRhs the local right-hand side vector
+   */
+  template< typename POLICY >
+  static void
+  createAndLaunch( integer const numComps,
+                   integer const numPhases,
+                   globalIndex const rankOffset,
+                   integer const useTotalMassEquation,
+                   integer const useSimpleAccumulation,
+                   string const dofKey,
+                   ElementSubRegionBase const & subRegion,
+                   constitutive::MultiFluidBase const & fluid,
+                   constitutive::CoupledSolidBase const & solid,
+                   CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                   arrayView1d< real64 > const & localRhs )
+  {
+    internal::kernelLaunchSelectorCompSwitch( numComps, [&] ( auto NC )
+    {
+      integer constexpr NUM_COMP = NC();
+      integer constexpr NUM_DOF = NC()+1;
+
+      BitFlags< KernelFlags > KernelFlags;
+      if( useTotalMassEquation )
+        KernelFlags.set( KernelFlags::TotalMassEquation );
+      if( useSimpleAccumulation )
+        KernelFlags.set( KernelFlags::SimpleAccumulation );
+
+      AccumulationZFormulationKernel< NUM_COMP, NUM_DOF > kernel( numPhases, rankOffset, dofKey, subRegion,
+                                                                  fluid, solid, localMatrix, localRhs, KernelFlags );
+      AccumulationZFormulationKernel< NUM_COMP, NUM_DOF >::template launch< POLICY >( subRegion.size(), kernel );
+    } );
+  }
+
+};
+
+} // namespace isothermalCompositionalMultiphaseBaseKernels
+
+} // namespace geos
+
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_ACCUMULATIONZFORMULATIONKERNEL_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/DirichletFluxComputeZFormulationKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/DirichletFluxComputeZFormulationKernel.hpp
new file mode 100644
index 0000000000..6e0455f4cc
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/DirichletFluxComputeZFormulationKernel.hpp
@@ -0,0 +1,544 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file DirichletFluxComputeZFormulationKernel.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_DIRICHLETFLUXCOMPUTEZFORMULATIONKERNEL_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_DIRICHLETFLUXCOMPUTEZFORMULATIONKERNEL_HPP
+
+#include "codingUtilities/Utilities.hpp"
+#include "common/DataLayouts.hpp"
+#include "common/DataTypes.hpp"
+#include "common/GEOS_RAJA_Interface.hpp"
+#include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
+#include "constitutive/fluid/multifluid/MultiFluidSelector.hpp"
+#include "finiteVolume/BoundaryStencil.hpp"
+#include "mesh/ElementRegionManager.hpp"
+#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
+#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp"
+#include "physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp"
+#include "physicsSolvers/fluidFlow/StencilAccessors.hpp"
+
+
+namespace geos
+{
+
+namespace isothermalCompositionalMultiphaseFVMKernels
+{
+
+/******************************** DirichletFluxComputeZFormulationKernel ********************************/
+
+/**
+ * @class DirichletFluxComputeZFormulationKernel
+ * @tparam NUM_COMP number of fluid components
+ * @tparam NUM_DOF number of degrees of freedom
+ * @tparam FLUIDWRAPPER the type of the fluid wrapper
+ * @brief Define the interface for the assembly kernel in charge of Dirichlet face flux terms
+ */
+template< integer NUM_COMP, integer NUM_DOF, typename FLUIDWRAPPER >
+class DirichletFluxComputeZFormulationKernel : public FluxComputeKernel< NUM_COMP,
+                                                                         NUM_DOF,
+                                                                         BoundaryStencilWrapper >
+{
+public:
+
+  /**
+   * @brief The type for element-based data. Consists entirely of ArrayView's.
+   *
+   * Can be converted from ElementRegionManager::ElementViewConstAccessor
+   * by calling .toView() or .toViewConst() on an accessor instance
+   */
+  template< typename VIEWTYPE >
+  using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
+
+  using AbstractBase = isothermalCompositionalMultiphaseFVMKernels::FluxComputeKernelBase;
+  using DofNumberAccessor = AbstractBase::DofNumberAccessor;
+  using CompFlowAccessors = AbstractBase::CompFlowAccessors;
+  using MultiFluidAccessors = AbstractBase::MultiFluidAccessors;
+  using CapPressureAccessors = AbstractBase::CapPressureAccessors;
+  using PermeabilityAccessors = AbstractBase::PermeabilityAccessors;
+
+  using AbstractBase::m_dt;
+  using AbstractBase::m_numPhases;
+  using AbstractBase::m_rankOffset;
+  using AbstractBase::m_dofNumber;
+  using AbstractBase::m_ghostRank;
+  using AbstractBase::m_gravCoef;
+  using AbstractBase::m_pres;
+  using AbstractBase::m_phaseCompFrac;
+  using AbstractBase::m_dPhaseCompFrac;
+  using AbstractBase::m_localMatrix;
+  using AbstractBase::m_localRhs;
+  using AbstractBase::m_kernelFlags;
+
+  using Base = isothermalCompositionalMultiphaseFVMKernels::FluxComputeKernel< NUM_COMP, NUM_DOF, BoundaryStencilWrapper >;
+  using Base::numComp;
+  using Base::numDof;
+  using Base::numEqn;
+  using Base::m_stencilWrapper;
+  using Base::m_phaseMob;
+  using Base::m_dPhaseMob;
+  using Base::m_phaseMassDens;
+  using Base::m_dPhaseMassDens;
+  using Base::m_permeability;
+  using Base::m_dPerm_dPres;
+  using Base::m_seri;
+  using Base::m_sesri;
+  using Base::m_sei;
+
+  /**
+   * @brief Constructor for the kernel interface
+   * @param[in] numPhases the number of fluid phases
+   * @param[in] rankOffset the offset of my MPI rank
+   * @param[in] faceManager the face manager
+   * @param[in] stencilWrapper reference to the stencil wrapper
+   * @param[in] fluidWrapper reference to the fluid wrapper
+   * @param[in] dofNumberAccessor
+   * @param[in] compFlowAccessors
+   * @param[in] multiFluidAccessors
+   * @param[in] capPressureAccessors
+   * @param[in] permeabilityAccessors
+   * @param[in] dt time step size
+   * @param[inout] localMatrix the local CRS matrix
+   * @param[inout] localRhs the local right-hand side vector
+   * @param[in] kernelFlags flags packed together
+   */
+  DirichletFluxComputeZFormulationKernel( integer const numPhases,
+                                          globalIndex const rankOffset,
+                                          FaceManager const & faceManager,
+                                          BoundaryStencilWrapper const & stencilWrapper,
+                                          FLUIDWRAPPER const & fluidWrapper,
+                                          DofNumberAccessor const & dofNumberAccessor,
+                                          CompFlowAccessors const & compFlowAccessors,
+                                          MultiFluidAccessors const & multiFluidAccessors,
+                                          CapPressureAccessors const & capPressureAccessors,
+                                          PermeabilityAccessors const & permeabilityAccessors,
+                                          real64 const dt,
+                                          CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                                          arrayView1d< real64 > const & localRhs,
+                                          BitFlags< KernelFlags > kernelFlags )
+    : Base( numPhases,
+            rankOffset,
+            stencilWrapper,
+            dofNumberAccessor,
+            compFlowAccessors,
+            multiFluidAccessors,
+            capPressureAccessors,
+            permeabilityAccessors,
+            dt,
+            localMatrix,
+            localRhs,
+            kernelFlags ),
+    m_facePres( faceManager.getField< fields::flow::facePressure >() ),
+    m_faceTemp( faceManager.getField< fields::flow::faceTemperature >() ),
+    m_faceCompFrac( faceManager.getField< fields::flow::faceGlobalCompFraction >() ),
+    m_faceGravCoef( faceManager.getField< fields::flow::gravityCoefficient >() ),
+    m_fluidWrapper( fluidWrapper )
+  {}
+
+  /**
+   * @struct StackVariables
+   * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack
+   */
+  struct StackVariables
+  {
+public:
+
+    /**
+     * @brief Constructor for the stack variables
+     * @param[in] size size of the stencil for this connection
+     * @param[in] numElems number of elements for this connection
+     */
+    GEOS_HOST_DEVICE
+    StackVariables( localIndex const GEOS_UNUSED_PARAM( size ),
+                    localIndex GEOS_UNUSED_PARAM( numElems )) {}
+
+    // Transmissibility
+    real64 transmissibility = 0.0;
+
+    // Component fluxes and derivatives
+
+    /// Component fluxes
+    real64 compFlux[numComp]{};
+    /// Derivatives of component fluxes wrt pressure
+    real64 dCompFlux_dP[numComp]{};
+    /// Derivatives of component fluxes wrt component densities
+    real64 dCompFlux_dC[numComp][numComp]{};
+
+    // Local degrees of freedom and local residual/jacobian
+
+    /// Indices of the matrix rows/columns corresponding to the dofs in this face
+    globalIndex dofColIndices[numDof]{};
+
+    /// Storage for the face local residual vector
+    real64 localFlux[numEqn]{};
+    /// Storage for the face local Jacobian matrix
+    real64 localFluxJacobian[numEqn][numDof]{};
+
+  };
+
+
+  /**
+   * @brief Performs the setup phase for the kernel.
+   * @param[in] iconn the connection index
+   * @param[in] stack the stack variables
+   */
+  GEOS_HOST_DEVICE
+  void setup( localIndex const iconn,
+              StackVariables & stack ) const
+  {
+    globalIndex const offset =
+      m_dofNumber[m_seri( iconn, BoundaryStencil::Order::ELEM )][m_sesri( iconn, BoundaryStencil::Order::ELEM )][m_sei( iconn, BoundaryStencil::Order::ELEM )];
+
+    for( integer jdof = 0; jdof < numDof; ++jdof )
+    {
+      stack.dofColIndices[jdof] = offset + jdof;
+    }
+  }
+
+
+  /**
+   * @brief Compute the local Dirichlet face flux contributions to the residual and Jacobian
+   * @tparam FUNC the type of the function that can be used to customize the computation of the phase fluxes
+   * @param[in] iconn the connection index
+   * @param[inout] stack the stack variables
+   * @param[in] compFluxKernelOp the function used to customize the computation of the component fluxes
+   */
+  template< typename FUNC = NoOpFunc >
+  GEOS_HOST_DEVICE
+  void computeFlux( localIndex const iconn,
+                    StackVariables & stack,
+                    FUNC && compFluxKernelOp = NoOpFunc{} ) const
+  {
+    using Deriv = constitutive::multifluid::DerivativeOffset;
+    using Order = BoundaryStencil::Order;
+
+    localIndex const er  = m_seri( iconn, Order::ELEM );
+    localIndex const esr = m_sesri( iconn, Order::ELEM );
+    localIndex const ei  = m_sei( iconn, Order::ELEM );
+    localIndex const kf  = m_sei( iconn, Order::FACE );
+
+    // Step 1: compute the transmissibility at the boundary face
+
+    real64 dTrans_dPerm[3]{};
+    m_stencilWrapper.computeWeights( iconn,
+                                     m_permeability,
+                                     stack.transmissibility,
+                                     dTrans_dPerm );
+    real64 const dTrans_dPres = LvArray::tensorOps::AiBi< 3 >( dTrans_dPerm, m_dPerm_dPres[er][esr][ei][0] );
+
+    // Step 2: compute the fluid properties on the face
+    // This is needed to get the phase mass density and the phase comp fraction at the face
+    // Because we approximate the face mobility using the total element mobility
+
+    StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseFrac( 1, 1, m_numPhases );
+    StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseDens( 1, 1, m_numPhases );
+    StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseMassDens( 1, 1, m_numPhases );
+    StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseVisc( 1, 1, m_numPhases );
+    StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseEnthalpy( 1, 1, m_numPhases );
+    StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseInternalEnergy( 1, 1, m_numPhases );
+    StackArray< real64, 4, constitutive::MultiFluidBase::MAX_NUM_PHASES * NUM_COMP,
+                constitutive::multifluid::LAYOUT_PHASE_COMP > facePhaseCompFrac( 1, 1, m_numPhases, NUM_COMP );
+    real64 faceTotalDens = 0.0;
+
+    constitutive::MultiFluidBase::KernelWrapper::computeValues( m_fluidWrapper,
+                                                                m_facePres[kf],
+                                                                m_faceTemp[kf],
+                                                                m_faceCompFrac[kf],
+                                                                facePhaseFrac[0][0],
+                                                                facePhaseDens[0][0],
+                                                                facePhaseMassDens[0][0],
+                                                                facePhaseVisc[0][0],
+                                                                facePhaseEnthalpy[0][0],
+                                                                facePhaseInternalEnergy[0][0],
+                                                                facePhaseCompFrac[0][0],
+                                                                faceTotalDens );
+
+    // Step 3: loop over phases, compute and upwind phase flux and sum contributions to each component's flux
+
+    for( integer ip = 0; ip < m_numPhases; ++ip )
+    {
+
+      // working variables
+      real64 dDensMean_dC[numComp]{};
+      real64 dF_dC[numComp]{};
+
+      real64 phaseFlux = 0.0; // for the lambda
+      real64 dPhaseFlux_dP = 0.0;
+      real64 dPhaseFlux_dC[numComp]{};
+
+
+      // Step 3.1: compute the average phase mass density at the face
+      // average density and derivatives
+      real64 const densMean = 0.5 * ( m_phaseMassDens[er][esr][ei][0][ip] + facePhaseMassDens[0][0][ip] );
+      real64 const dDensMean_dP = 0.5 * m_dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP];
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        dDensMean_dC[jc] = 0.5 * m_dPhaseMassDens[er][esr][ei][0][ip][Deriv::dC+jc];
+      }
+
+
+      // Step 3.2: compute the (TPFA) potential difference at the face
+
+      real64 const gravTimesDz = m_gravCoef[er][esr][ei] - m_faceGravCoef[kf];
+      real64 const potDif = m_pres[er][esr][ei] - m_facePres[kf] - densMean * gravTimesDz;
+      real64 const f = stack.transmissibility * potDif;
+      real64 const dF_dP = stack.transmissibility * ( 1.0 - dDensMean_dP * gravTimesDz ) + dTrans_dPres * potDif;
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        dF_dC[jc] = -stack.transmissibility * dDensMean_dC[jc] * gravTimesDz;
+      }
+
+      // Step 3.3: computation of the mobility
+      // We do that before the if/else statement to be able to pass it to the compFluxOpKernel
+
+      // recomputing the exact mobility at the face would be quite complex, as it would require:
+      //   1) computing the saturation
+      //   2) computing the relperm
+      //   3) computing the mobility as \lambda_p = \rho_p kr_p( S_p ) / \mu_p
+      // the second step in particular would require yet another dispatch to get the relperm model
+      // so, for simplicity, we approximate the face mobility as
+      //    \lambda^approx_p = \rho_p S_p / \mu_p
+      //                     = \rho_p ( (nu_p / rho_p) * rho_t ) / \mu_p (plugging the expression of saturation)
+      //                     = \nu_p * rho_t / \mu_p
+      // fortunately, we don't need the derivatives
+      real64 const facePhaseMob = ( facePhaseFrac[0][0][ip] > 0.0 )
+  ? facePhaseFrac[0][0][ip] * faceTotalDens / facePhaseVisc[0][0][ip]
+  : 0.0;
+
+      // *** upwinding ***
+      // Step 3.4: upwinding based on the sign of the phase potential gradient
+      // It is easier to hard-code the if/else because it is difficult to address elem and face variables in a uniform way
+
+      if( potDif >= 0 ) // the element is upstream
+      {
+
+        // compute the phase flux and derivatives using the element mobility
+        phaseFlux = m_phaseMob[er][esr][ei][ip] * f;
+        dPhaseFlux_dP = m_phaseMob[er][esr][ei][ip] * dF_dP + m_dPhaseMob[er][esr][ei][ip][Deriv::dP] * f;
+        for( integer jc = 0; jc < numComp; ++jc )
+        {
+          dPhaseFlux_dC[jc] =
+            m_phaseMob[er][esr][ei][ip] * dF_dC[jc] + m_dPhaseMob[er][esr][ei][ip][Deriv::dC+jc] * f;
+        }
+
+        // slice some constitutive arrays to avoid too much indexing in component loop
+        arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE_COMP-3 > phaseCompFracSub =
+          m_phaseCompFrac[er][esr][ei][0][ip];
+        arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC-3 > dPhaseCompFracSub =
+          m_dPhaseCompFrac[er][esr][ei][0][ip];
+
+        // compute component fluxes and derivatives using element composition
+        for( integer ic = 0; ic < numComp; ++ic )
+        {
+          real64 const ycp = phaseCompFracSub[ic];
+          stack.compFlux[ic] += phaseFlux * ycp;
+          stack.dCompFlux_dP[ic] += dPhaseFlux_dP * ycp + phaseFlux * dPhaseCompFracSub[ic][Deriv::dP];
+
+          for( integer jc = 0; jc < numComp; ++jc )
+          {
+            stack.dCompFlux_dC[ic][jc] += dPhaseFlux_dC[jc] * ycp + phaseFlux * dPhaseCompFracSub[ic][Deriv::dC+jc];
+          }
+        }
+
+      }
+      else // the face is upstream
+      {
+
+        // compute the phase flux and derivatives using the approximated face mobility
+        // we only have to take derivatives of the potential gradient in this case
+        phaseFlux = facePhaseMob * f;
+        dPhaseFlux_dP = facePhaseMob * dF_dP;
+        for( integer jc = 0; jc < numComp; ++jc )
+        {
+          dPhaseFlux_dC[jc] = facePhaseMob * dF_dC[jc];
+        }
+
+        // compute component fluxes and derivatives using the face composition
+        for( integer ic = 0; ic < numComp; ++ic )
+        {
+          real64 const ycp = facePhaseCompFrac[0][0][ip][ic];
+          stack.compFlux[ic] += phaseFlux * ycp;
+          stack.dCompFlux_dP[ic] += dPhaseFlux_dP * ycp;
+          for( integer jc = 0; jc < numComp; ++jc )
+          {
+            stack.dCompFlux_dC[ic][jc] += dPhaseFlux_dC[jc] * ycp;
+          }
+        }
+      }
+
+      // call the lambda in the phase loop to allow the reuse of the phase fluxes and their derivatives
+      // possible use: assemble the derivatives wrt temperature, and the flux term of the energy equation for this phase
+      compFluxKernelOp( ip, er, esr, ei, kf, f,
+                        facePhaseMob, facePhaseEnthalpy[0][0], facePhaseCompFrac[0][0],
+                        phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
+
+    }
+
+    // *** end of upwinding
+
+    // Step 4: populate local flux vector and derivatives
+    for( integer ic = 0; ic < numComp; ++ic )
+    {
+      stack.localFlux[ic]            = m_dt * stack.compFlux[ic];
+      stack.localFluxJacobian[ic][0] = m_dt * stack.dCompFlux_dP[ic];
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        stack.localFluxJacobian[ic][jc+1] = m_dt * stack.dCompFlux_dC[ic][jc];
+      }
+    }
+  }
+
+  /**
+   * @brief Performs the complete phase for the kernel.
+   * @param[in] iconn the connection index
+   * @param[inout] stack the stack variables
+   */
+  template< typename FUNC = NoOpFunc >
+  GEOS_HOST_DEVICE
+  void complete( localIndex const iconn,
+                 StackVariables & stack,
+                 FUNC && assemblyKernelOp = NoOpFunc{} ) const
+  {
+    using namespace compositionalMultiphaseUtilities;
+    using Order = BoundaryStencil::Order;
+
+    if( AbstractBase::m_kernelFlags.isSet( KernelFlags::TotalMassEquation ) )
+    {
+      // Apply equation/variable change transformation(s)
+      real64 work[numDof]{};
+      shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numDof, stack.localFluxJacobian, work );
+      shiftElementsAheadByOneAndReplaceFirstElementWithSum( numComp, stack.localFlux );
+    }
+
+    // add contribution to residual and jacobian into:
+    // - the component mass balance equations (i = 0 to i = numComp-1)
+    // note that numDof includes derivatives wrt temperature if this class is derived in ThermalKernels
+    if( m_ghostRank[m_seri( iconn, Order::ELEM )][m_sesri( iconn, Order::ELEM )][m_sei( iconn, Order::ELEM )] < 0 )
+    {
+      globalIndex const globalRow = m_dofNumber[m_seri( iconn, Order::ELEM )][m_sesri( iconn, Order::ELEM )][m_sei( iconn, Order::ELEM )];
+      localIndex const localRow = LvArray::integerConversion< localIndex >( globalRow - m_rankOffset );
+      GEOS_ASSERT_GE( localRow, 0 );
+      GEOS_ASSERT_GT( AbstractBase::m_localMatrix.numRows(), localRow + numComp );
+
+      for( integer ic = 0; ic < numComp; ++ic )
+      {
+        RAJA::atomicAdd( parallelDeviceAtomic{}, &AbstractBase::m_localRhs[localRow + ic], stack.localFlux[ic] );
+        AbstractBase::m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >
+          ( localRow + ic,
+          stack.dofColIndices,
+          stack.localFluxJacobian[ic],
+          numDof );
+      }
+
+      // call the lambda to assemble additional terms, such as thermal terms
+      assemblyKernelOp( localRow );
+    }
+  }
+
+protected:
+
+  /// Views on face pressure, temperature, and composition
+  arrayView1d< real64 const > const m_facePres;
+  arrayView1d< real64 const > const m_faceTemp;
+  arrayView2d< real64 const, compflow::USD_COMP > const m_faceCompFrac;
+
+  /// View on the face gravity coefficient
+  arrayView1d< real64 const > const m_faceGravCoef;
+
+  /// Reference to the fluid wrapper
+  FLUIDWRAPPER const m_fluidWrapper;
+
+};
+
+
+/**
+ * @class DirichletFluxComputeZFormulationKernelFactory
+ */
+class DirichletFluxComputeZFormulationKernelFactory
+{
+public:
+
+  /**
+   * @brief Create a new kernel and launch
+   * @tparam POLICY the policy used in the RAJA kernel
+   * @param[in] numComps the number of fluid components
+   * @param[in] numPhases the number of fluid phases
+   * @param[in] rankOffset the offset of my MPI rank
+   * @param[in] dofKey string to get the element degrees of freedom numbers
+   * @param[in] solverName name of the solver (to name accessors)
+   * @param[in] faceManager reference to the face manager
+   * @param[in] elemManager reference to the element region manager
+   * @param[in] stencilWrapper reference to the boundary stencil wrapper
+   * @param[in] fluidBase the multifluid constitutive model
+   * @param[in] dt time step size
+   * @param[inout] localMatrix the local CRS matrix
+   * @param[inout] localRhs the local right-hand side vector
+   */
+  template< typename POLICY >
+  static void
+  createAndLaunch( integer const numComps,
+                   integer const numPhases,
+                   globalIndex const rankOffset,
+                   integer const useTotalMassEquation,
+                   string const & dofKey,
+                   string const & solverName,
+                   FaceManager const & faceManager,
+                   ElementRegionManager const & elemManager,
+                   BoundaryStencilWrapper const & stencilWrapper,
+                   constitutive::MultiFluidBase & fluidBase,
+                   real64 const dt,
+                   CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                   arrayView1d< real64 > const & localRhs )
+  {
+    constitutive::constitutiveComponentUpdatePassThru( fluidBase, numComps, [&]( auto & fluid, auto NC )
+    {
+      using FluidType = TYPEOFREF( fluid );
+      typename FluidType::KernelWrapper const fluidWrapper = fluid.createKernelWrapper();
+
+      integer constexpr NUM_COMP = NC();
+      integer constexpr NUM_DOF = NC() + 1;
+
+      ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > dofNumberAccessor =
+        elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey );
+      dofNumberAccessor.setName( solverName + "/accessors/" + dofKey );
+
+      // for now, we neglect capillary pressure in the kernel
+      BitFlags< KernelFlags > kernelFlags;
+      if( useTotalMassEquation )
+        kernelFlags.set( KernelFlags::TotalMassEquation );
+
+      using kernelType = DirichletFluxComputeZFormulationKernel< NUM_COMP, NUM_DOF, typename FluidType::KernelWrapper >;
+      typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName );
+      typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName );
+      typename kernelType::CapPressureAccessors capPressureAccessors( elemManager, solverName );
+      typename kernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName );
+
+      kernelType kernel( numPhases, rankOffset, faceManager, stencilWrapper, fluidWrapper,
+                         dofNumberAccessor, compFlowAccessors, multiFluidAccessors, capPressureAccessors, permeabilityAccessors,
+                         dt, localMatrix, localRhs, kernelFlags );
+      kernelType::template launch< POLICY >( stencilWrapper.size(), kernel );
+    } );
+  }
+};
+
+} // namespace isothermalCompositionalMultiphaseFVMKernels
+
+} // namespace geos
+
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_DIRICHLETFLUXCOMPUTEZFORMULATIONKERNEL_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/FluxComputeZFormulationKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/FluxComputeZFormulationKernel.hpp
new file mode 100644
index 0000000000..0286d1e674
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/FluxComputeZFormulationKernel.hpp
@@ -0,0 +1,528 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file FluxComputeZFormulationKernel.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_FLUXCOMPUTEZFORMULATIONKERNEL_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_FLUXCOMPUTEZFORMULATIONKERNEL_HPP
+
+#include "physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp"
+
+#include "codingUtilities/Utilities.hpp"
+#include "common/DataLayouts.hpp"
+#include "common/DataTypes.hpp"
+#include "common/GEOS_RAJA_Interface.hpp"
+#include "constitutive/fluid/multifluid/MultiFluidFields.hpp"
+#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
+#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp"
+#include "physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp"
+#include "physicsSolvers/fluidFlow/StencilAccessors.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/KernelLaunchSelectors.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PPUPhaseFluxZFormulation.hpp"
+
+namespace geos
+{
+
+namespace isothermalCompositionalMultiphaseFVMKernels
+{
+
+/**
+ * @class FluxComputeZFormulationKernel
+ * @tparam NUM_COMP number of fluid components
+ * @tparam NUM_DOF number of degrees of freedom
+ * @tparam STENCILWRAPPER the type of the stencil wrapper
+ * @brief Define the interface for the assembly kernel in charge of flux terms
+ */
+template< integer NUM_COMP, integer NUM_DOF, typename STENCILWRAPPER >
+class FluxComputeZFormulationKernel : public FluxComputeKernelBase
+{
+public:
+
+  /// Compile time value for the number of components
+  static constexpr integer numComp = NUM_COMP;
+
+  /// Compute time value for the number of degrees of freedom
+  static constexpr integer numDof = NUM_DOF;
+
+  /// Compute time value for the number of equations (all of them, except the volume balance equation)
+  static constexpr integer numEqn = NUM_DOF-1;
+
+  /// Maximum number of elements at the face
+  static constexpr localIndex maxNumElems = STENCILWRAPPER::maxNumPointsInFlux;
+
+  /// Maximum number of connections at the face
+  static constexpr localIndex maxNumConns = STENCILWRAPPER::maxNumConnections;
+
+  /// Maximum number of points in the stencil
+  static constexpr localIndex maxStencilSize = STENCILWRAPPER::maxStencilSize;
+
+  /// Number of flux support points (hard-coded for TFPA)
+  static constexpr integer numFluxSupportPoints = 2;
+
+  /**
+   * @brief Constructor for the kernel interface
+   * @param[in] numPhases the number of fluid phases
+   * @param[in] rankOffset the offset of my MPI rank
+   * @param[in] stencilWrapper reference to the stencil wrapper
+   * @param[in] dofNumberAccessor
+   * @param[in] compFlowAccessors
+   * @param[in] multiFluidAccessors
+   * @param[in] capPressureAccessors
+   * @param[in] permeabilityAccessors
+   * @param[in] dt time step size
+   * @param[inout] localMatrix the local CRS matrix
+   * @param[inout] localRhs the local right-hand side vector
+   * @param[in] kernelFlags flags packed together
+   */
+  FluxComputeZFormulationKernel( integer const numPhases,
+                                 globalIndex const rankOffset,
+                                 STENCILWRAPPER const & stencilWrapper,
+                                 DofNumberAccessor const & dofNumberAccessor,
+                                 CompFlowAccessors const & compFlowAccessors,
+                                 MultiFluidAccessors const & multiFluidAccessors,
+                                 CapPressureAccessors const & capPressureAccessors,
+                                 PermeabilityAccessors const & permeabilityAccessors,
+                                 real64 const dt,
+                                 CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                                 arrayView1d< real64 > const & localRhs,
+                                 BitFlags< KernelFlags > kernelFlags )
+    : FluxComputeKernelBase( numPhases,
+                             rankOffset,
+                             dofNumberAccessor,
+                             compFlowAccessors,
+                             multiFluidAccessors,
+                             dt,
+                             localMatrix,
+                             localRhs,
+                             kernelFlags ),
+    m_permeability( permeabilityAccessors.get( fields::permeability::permeability {} ) ),
+    m_dPerm_dPres( permeabilityAccessors.get( fields::permeability::dPerm_dPressure {} ) ),
+    m_phaseMob( compFlowAccessors.get( fields::flow::phaseMobility {} ) ),
+    m_dPhaseMob( compFlowAccessors.get( fields::flow::dPhaseMobility {} ) ),
+    m_phaseMassDens( multiFluidAccessors.get( fields::multifluid::phaseMassDensity {} ) ),
+    m_dPhaseMassDens( multiFluidAccessors.get( fields::multifluid::dPhaseMassDensity {} ) ),
+    m_phaseCapPressure( capPressureAccessors.get( fields::cappres::phaseCapPressure {} ) ),
+    m_dPhaseCapPressure_dPhaseVolFrac( capPressureAccessors.get( fields::cappres::dPhaseCapPressure_dPhaseVolFraction {} ) ),
+    m_stencilWrapper( stencilWrapper ),
+    m_seri( stencilWrapper.getElementRegionIndices() ),
+    m_sesri( stencilWrapper.getElementSubRegionIndices() ),
+    m_sei( stencilWrapper.getElementIndices() )
+  { }
+
+  /**
+   * @struct StackVariables
+   * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack
+   */
+  struct StackVariables
+  {
+public:
+
+    /**
+     * @brief Constructor for the stack variables
+     * @param[in] size size of the stencil for this connection
+     * @param[in] numElems number of elements for this connection
+     */
+    GEOS_HOST_DEVICE
+    StackVariables( localIndex const size, localIndex numElems )
+      : stencilSize( size ),
+      numConnectedElems( numElems ),
+      dofColIndices( size * numDof ),
+      localFlux( numElems * numEqn ),
+      localFluxJacobian( numElems * numEqn, size * numDof )
+    {}
+
+    // Stencil information
+
+    /// Stencil size for a given connection
+    localIndex const stencilSize;
+    /// Number of elements connected at a given connection
+    localIndex const numConnectedElems;
+
+    // Transmissibility and derivatives
+
+    /// Transmissibility
+    real64 transmissibility[maxNumConns][numFluxSupportPoints]{};
+    /// Derivatives of transmissibility with respect to pressure
+    real64 dTrans_dPres[maxNumConns][numFluxSupportPoints]{};
+
+    // Local degrees of freedom and local residual/jacobian
+
+    /// Indices of the matrix rows/columns corresponding to the dofs in this face
+    stackArray1d< globalIndex, maxNumElems * numDof > dofColIndices;
+
+    /// Storage for the face local residual vector (all equations except volume balance)
+    stackArray1d< real64, maxNumElems * numEqn > localFlux;
+    /// Storage for the face local Jacobian matrix
+    stackArray2d< real64, maxNumElems * numEqn * maxStencilSize * numDof > localFluxJacobian;
+  };
+
+
+  /**
+   * @brief Getter for the stencil size at this connection
+   * @param[in] iconn the connection index
+   * @return the size of the stencil at this connection
+   */
+  GEOS_HOST_DEVICE
+  inline
+  localIndex stencilSize( localIndex const iconn ) const { return m_sei[iconn].size(); }
+
+  /**
+   * @brief Getter for the number of elements at this connection
+   * @param[in] iconn the connection index
+   * @return the number of elements at this connection
+   */
+  GEOS_HOST_DEVICE
+  inline
+  localIndex numPointsInFlux( localIndex const iconn ) const { return m_stencilWrapper.numPointsInFlux( iconn ); }
+
+
+  /**
+   * @brief Performs the setup phase for the kernel.
+   * @param[in] iconn the connection index
+   * @param[in] stack the stack variables
+   */
+  GEOS_HOST_DEVICE
+  inline
+  void setup( localIndex const iconn,
+              StackVariables & stack ) const
+  {
+    // set degrees of freedom indices for this face
+    for( integer i = 0; i < stack.stencilSize; ++i )
+    {
+      globalIndex const offset = m_dofNumber[m_seri( iconn, i )][m_sesri( iconn, i )][m_sei( iconn, i )];
+
+      for( integer jdof = 0; jdof < numDof; ++jdof )
+      {
+        stack.dofColIndices[i * numDof + jdof] = offset + jdof;
+      }
+    }
+  }
+
+  /**
+   * @brief Compute the local flux contributions to the residual and Jacobian
+   * @tparam FUNC the type of the function that can be used to customize the computation of the phase fluxes
+   * @param[in] iconn the connection index
+   * @param[inout] stack the stack variables
+   * @param[in] compFluxKernelOp the function used to customize the computation of the component fluxes
+   */
+  template< typename FUNC = NoOpFunc >
+  GEOS_HOST_DEVICE
+  inline
+  void computeFlux( localIndex const iconn,
+                    StackVariables & stack,
+                    FUNC && compFluxKernelOp = NoOpFunc{} ) const
+  {
+
+    // first, compute the transmissibilities at this face
+    m_stencilWrapper.computeWeights( iconn,
+                                     m_permeability,
+                                     m_dPerm_dPres,
+                                     stack.transmissibility,
+                                     stack.dTrans_dPres );
+
+
+    localIndex k[numFluxSupportPoints];
+    localIndex connectionIndex = 0;
+    for( k[0] = 0; k[0] < stack.numConnectedElems; ++k[0] )
+    {
+      for( k[1] = k[0] + 1; k[1] < stack.numConnectedElems; ++k[1] )
+      {
+        /// cell indices
+        localIndex const seri[numFluxSupportPoints]  = {m_seri( iconn, k[0] ), m_seri( iconn, k[1] )};
+        localIndex const sesri[numFluxSupportPoints] = {m_sesri( iconn, k[0] ), m_sesri( iconn, k[1] )};
+        localIndex const sei[numFluxSupportPoints]   = {m_sei( iconn, k[0] ), m_sei( iconn, k[1] )};
+
+        // clear working arrays
+        real64 compFlux[numComp]{};
+        real64 dCompFlux_dP[numFluxSupportPoints][numComp]{};
+        real64 dCompFlux_dC[numFluxSupportPoints][numComp][numComp]{};
+
+        real64 const trans[numFluxSupportPoints] = { stack.transmissibility[connectionIndex][0],
+                                                     stack.transmissibility[connectionIndex][1] };
+
+        real64 const dTrans_dPres[numFluxSupportPoints] = { stack.dTrans_dPres[connectionIndex][0],
+                                                            stack.dTrans_dPres[connectionIndex][1] };
+
+        //***** calculation of flux *****
+        // loop over phases, compute and upwind phase flux and sum contributions to each component's flux
+        for( integer ip = 0; ip < m_numPhases; ++ip )
+        {
+          // create local work arrays
+          real64 potGrad = 0.0;
+          real64 phaseFlux = 0.0;
+          real64 dPhaseFlux_dP[numFluxSupportPoints]{};
+          real64 dPhaseFlux_dC[numFluxSupportPoints][numComp]{};
+
+          localIndex k_up = -1;
+
+          isothermalCompositionalMultiphaseFVMKernelUtilities::PPUPhaseFluxZFormulation::compute< numComp, numFluxSupportPoints >
+            ( m_numPhases,
+            ip,
+            m_kernelFlags.isSet( KernelFlags::CapPressure ),
+            m_kernelFlags.isSet( KernelFlags::NewGravity ),
+            seri, sesri, sei,
+            trans,
+            dTrans_dPres,
+            m_pres,
+            m_gravCoef,
+            m_phaseMob, m_dPhaseMob,
+            m_phaseVolFrac, m_dPhaseVolFrac,
+            m_phaseCompFrac, m_dPhaseCompFrac,
+            m_phaseMassDens, m_dPhaseMassDens,
+            m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac,
+            k_up,
+            potGrad,
+            phaseFlux,
+            dPhaseFlux_dP,
+            dPhaseFlux_dC,
+            compFlux,
+            dCompFlux_dP,
+            dCompFlux_dC );
+
+          // call the lambda in the phase loop to allow the reuse of the phase fluxes and their derivatives
+          // possible use: assemble the derivatives wrt temperature, and the flux term of the energy equation for this phase
+          compFluxKernelOp( ip, k, seri, sesri, sei, connectionIndex,
+                            k_up, seri[k_up], sesri[k_up], sei[k_up], potGrad,
+                            phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
+
+        }                                 // loop over phases
+
+        /// populate local flux vector and derivatives
+        for( integer ic = 0; ic < numComp; ++ic )
+        {
+          integer const eqIndex0 = k[0] * numEqn + ic;
+          integer const eqIndex1 = k[1] * numEqn + ic;
+
+          stack.localFlux[eqIndex0]  +=  m_dt * compFlux[ic];
+          stack.localFlux[eqIndex1]  -=  m_dt * compFlux[ic];
+
+          for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+          {
+            localIndex const localDofIndexPres = k[ke] * numDof;
+            stack.localFluxJacobian[eqIndex0][localDofIndexPres] += m_dt * dCompFlux_dP[ke][ic];
+            stack.localFluxJacobian[eqIndex1][localDofIndexPres] -= m_dt * dCompFlux_dP[ke][ic];
+
+            for( integer jc = 0; jc < numComp; ++jc )
+            {
+              localIndex const localDofIndexComp = localDofIndexPres + jc + 1;
+              stack.localFluxJacobian[eqIndex0][localDofIndexComp] += m_dt * dCompFlux_dC[ke][ic][jc];
+              stack.localFluxJacobian[eqIndex1][localDofIndexComp] -= m_dt * dCompFlux_dC[ke][ic][jc];
+            }
+          }
+        }
+        connectionIndex++;
+      }   // loop over k[1]
+    }   // loop over k[0]
+
+  }
+
+  /**
+   * @brief Performs the complete phase for the kernel.
+   * @param[in] iconn the connection index
+   * @param[inout] stack the stack variables
+   */
+  template< typename FUNC = NoOpFunc >
+  GEOS_HOST_DEVICE
+  inline
+  void complete( localIndex const iconn,
+                 StackVariables & stack,
+                 FUNC && assemblyKernelOp = NoOpFunc{} ) const
+  {
+    using namespace compositionalMultiphaseUtilities;
+
+    if( m_kernelFlags.isSet( KernelFlags::TotalMassEquation ) )
+    {
+      // Apply equation/variable change transformation(s)
+      stackArray1d< real64, maxStencilSize * numDof > work( stack.stencilSize * numDof );
+      shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numEqn, numDof * stack.stencilSize, stack.numConnectedElems,
+                                                               stack.localFluxJacobian, work );
+      shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComp, numEqn, stack.numConnectedElems,
+                                                                 stack.localFlux );
+    }
+
+    // add contribution to residual and jacobian into:
+    // - the component mass balance equations (i = 0 to i = numComp-1)
+    // note that numDof includes derivatives wrt temperature if this class is derived in ThermalKernels
+    for( integer i = 0; i < stack.numConnectedElems; ++i )
+    {
+      if( m_ghostRank[m_seri( iconn, i )][m_sesri( iconn, i )][m_sei( iconn, i )] < 0 )
+      {
+        globalIndex const globalRow = m_dofNumber[m_seri( iconn, i )][m_sesri( iconn, i )][m_sei( iconn, i )];
+        localIndex const localRow = LvArray::integerConversion< localIndex >( globalRow - m_rankOffset );
+        GEOS_ASSERT_GE( localRow, 0 );
+        GEOS_ASSERT_GT( m_localMatrix.numRows(), localRow + numComp );
+
+        for( integer ic = 0; ic < numComp; ++ic )
+        {
+          RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[localRow + ic],
+                           stack.localFlux[i * numEqn + ic] );
+          m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >
+            ( localRow + ic,
+            stack.dofColIndices.data(),
+            stack.localFluxJacobian[i * numEqn + ic].dataIfContiguous(),
+            stack.stencilSize * numDof );
+        }
+
+        // call the lambda to assemble additional terms, such as thermal terms
+        assemblyKernelOp( i, localRow );
+      }
+    }
+  }
+
+  /**
+   * @brief Performs the kernel launch
+   * @tparam POLICY the policy used in the RAJA kernels
+   * @tparam KERNEL_TYPE the kernel type
+   * @param[in] numConnections the number of connections
+   * @param[inout] kernelComponent the kernel component providing access to setup/compute/complete functions and stack variables
+   */
+  template< typename POLICY, typename KERNEL_TYPE >
+  static void
+  launch( localIndex const numConnections,
+          KERNEL_TYPE const & kernelComponent )
+  {
+    GEOS_MARK_FUNCTION;
+    forAll< POLICY >( numConnections, [=] GEOS_HOST_DEVICE ( localIndex const iconn )
+    {
+      typename KERNEL_TYPE::StackVariables stack( kernelComponent.stencilSize( iconn ),
+                                                  kernelComponent.numPointsInFlux( iconn ) );
+
+      kernelComponent.setup( iconn, stack );
+      kernelComponent.computeFlux( iconn, stack );
+      kernelComponent.complete( iconn, stack );
+    } );
+  }
+
+protected:
+
+  /// Views on permeability
+  ElementViewConst< arrayView3d< real64 const > > const m_permeability;
+  ElementViewConst< arrayView3d< real64 const > > const m_dPerm_dPres;
+
+  /// Views on phase mobilities
+  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_phaseMob;
+  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const m_dPhaseMob;
+
+  /// Views on phase mass densities
+  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const m_phaseMassDens;
+  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const m_dPhaseMassDens;
+
+  /// Views on phase capillary pressure
+  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const m_phaseCapPressure;
+  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const m_dPhaseCapPressure_dPhaseVolFrac;
+
+  // Stencil information
+
+  /// Reference to the stencil wrapper
+  STENCILWRAPPER const m_stencilWrapper;
+
+  /// Connection to element maps
+  typename STENCILWRAPPER::IndexContainerViewConstType const m_seri;
+  typename STENCILWRAPPER::IndexContainerViewConstType const m_sesri;
+  typename STENCILWRAPPER::IndexContainerViewConstType const m_sei;
+
+};
+
+/**
+ * @class FluxComputeZFormulationKernelFactory
+ */
+class FluxComputeZFormulationKernelFactory
+{
+public:
+
+  /**
+   * @brief Create a new kernel and launch
+   * @tparam POLICY the policy used in the RAJA kernel
+   * @tparam STENCILWRAPPER the type of the stencil wrapper
+   * @param[in] numComps the number of fluid components
+   * @param[in] numPhases the number of fluid phases
+   * @param[in] rankOffset the offset of my MPI rank
+   * @param[in] dofKey string to get the element degrees of freedom numbers
+   * @param[in] hasCapPressure flag specifying whether capillary pressure is used or not
+   * @param[in] solverName name of the solver (to name accessors)
+   * @param[in] elemManager reference to the element region manager
+   * @param[in] stencilWrapper reference to the stencil wrapper
+   * @param[in] dt time step size
+   * @param[inout] localMatrix the local CRS matrix
+   * @param[inout] localRhs the local right-hand side vector
+   */
+  template< typename POLICY, typename STENCILWRAPPER >
+  static void
+  createAndLaunch( integer const numComps,
+                   integer const numPhases,
+                   globalIndex const rankOffset,
+                   string const & dofKey,
+                   integer const hasCapPressure,
+                   integer const useTotalMassEquation,
+                   integer const useNewGravity,
+                   UpwindingParameters upwindingParams,
+                   string const & solverName,
+                   ElementRegionManager const & elemManager,
+                   STENCILWRAPPER const & stencilWrapper,
+                   real64 const dt,
+                   CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                   arrayView1d< real64 > const & localRhs )
+  {
+    isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&]( auto NC )
+    {
+      integer constexpr NUM_COMP = NC();
+      integer constexpr NUM_DOF = NC() + 1;
+
+      ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > dofNumberAccessor =
+        elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey );
+      dofNumberAccessor.setName( solverName + "/accessors/" + dofKey );
+
+      BitFlags< KernelFlags > kernelFlags;
+      if( hasCapPressure )
+        kernelFlags.set( KernelFlags::CapPressure );
+      if( useTotalMassEquation )
+        kernelFlags.set( KernelFlags::TotalMassEquation );
+      if( useNewGravity )
+        kernelFlags.set( KernelFlags::NewGravity );
+      if( upwindingParams.upwindingScheme == UpwindingScheme::C1PPU &&
+          isothermalCompositionalMultiphaseFVMKernelUtilities::epsC1PPU > 0 )
+      {
+        GEOS_ERROR( "Z Formulation is currently not available for C1PPU" );
+      }
+      else if( upwindingParams.upwindingScheme == UpwindingScheme::IHU )
+      {
+        GEOS_ERROR( "Z Formulation is currently not available for IHU" );
+      }
+
+
+      using kernelType = FluxComputeZFormulationKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >;
+      typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName );
+      typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName );
+      typename kernelType::CapPressureAccessors capPressureAccessors( elemManager, solverName );
+      typename kernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName );
+
+      kernelType kernel( numPhases, rankOffset, stencilWrapper, dofNumberAccessor,
+                         compFlowAccessors, multiFluidAccessors, capPressureAccessors, permeabilityAccessors,
+                         dt, localMatrix, localRhs, kernelFlags );
+      kernelType::template launch< POLICY >( stencilWrapper.size(), kernel );
+    } );
+  }
+};
+
+} // namespace isothermalCompositionalMultiphaseFVMKernels
+
+} // namespace geos
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_FLUXCOMPUTEZFORMULATIONKERNEL_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PPUPhaseFluxZFormulation.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PPUPhaseFluxZFormulation.hpp
new file mode 100644
index 0000000000..cfb4ba9b12
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PPUPhaseFluxZFormulation.hpp
@@ -0,0 +1,161 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file PPUPhaseFluxZFormulation.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PPUPHASEFLUXZFORMULATION_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PPUPHASEFLUXZFORMULATION_HPP
+
+#include "common/DataLayouts.hpp"
+#include "common/DataTypes.hpp"
+#include "constitutive/fluid/multifluid/Layouts.hpp"
+#include "constitutive/capillaryPressure/layouts.hpp"
+#include "mesh/ElementRegionManager.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PotGradZFormulation.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseComponentFluxZFormulation.hpp"
+
+namespace geos
+{
+
+namespace isothermalCompositionalMultiphaseFVMKernelUtilities
+{
+
+template< typename VIEWTYPE >
+using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
+
+using Deriv = constitutive::multifluid::DerivativeOffset;
+
+struct PPUPhaseFluxZFormulation
+{
+  /**
+   * @brief Form the PhasePotentialUpwind from pressure gradient and gravitational head
+   * @tparam numComp number of components
+   * @tparam numFluxSupportPoints number of flux support points
+   * @param numPhase number of phases
+   * @param ip phase index
+   * @param hasCapPressure flag indicating if there is capillary pressure
+   * @param seri arraySlice of the stencil-implied element region index
+   * @param sesri arraySlice of the stencil-implied element subregion index
+   * @param sei arraySlice of the stencil-implied element index
+   * @param trans transmissibility at the connection
+   * @param dTrans_dPres derivative of transmissibility wrt pressure
+   * @param pres pressure
+   * @param gravCoef gravitational coefficient
+   * @param phaseMob phase mobility
+   * @param dPhaseMob derivative of phase mobility wrt pressure, temperature, comp density
+   * @param dPhaseVolFrac derivative of phase volume fraction wrt pressure, temperature, comp density
+   * @param phaseMassDens phase mass density
+   * @param dPhaseMassDens derivative of phase mass density wrt pressure, temperature, comp fraction
+   * @param phaseCapPressure phase capillary pressure
+   * @param dPhaseCapPressure_dPhaseVolFrac derivative of phase capillary pressure wrt phase volume fraction
+   * @param k_up uptream index for this phase
+   * @param potGrad potential gradient for this phase
+   * @param phaseFlux phase flux
+   * @param dPhaseFlux_dP derivative of phase flux wrt pressure
+   * @param dPhaseFlux_dC derivative of phase flux wrt comp density
+   */
+  template< integer numComp, integer numFluxSupportPoints >
+  GEOS_HOST_DEVICE
+  static void
+  compute( integer const numPhase,
+           integer const ip,
+           integer const hasCapPressure,
+           integer const useNewGravity,
+           localIndex const ( &seri )[numFluxSupportPoints],
+           localIndex const ( &sesri )[numFluxSupportPoints],
+           localIndex const ( &sei )[numFluxSupportPoints],
+           real64 const ( &trans )[2],
+           real64 const ( &dTrans_dPres )[2],
+           ElementViewConst< arrayView1d< real64 const > > const & pres,
+           ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
+           ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
+           ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
+           ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
+           ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
+           ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac,
+           ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac,
+           ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
+           ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
+           ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
+           ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
+           localIndex & k_up,
+           real64 & potGrad,
+           real64 ( &phaseFlux ),
+           real64 ( & dPhaseFlux_dP )[numFluxSupportPoints],
+           real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp],
+           real64 ( & compFlux )[numComp],
+           real64 ( & dCompFlux_dP )[numFluxSupportPoints][numComp],
+           real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp] )
+  {
+    real64 dPresGrad_dP[numFluxSupportPoints]{};
+    real64 dPresGrad_dC[numFluxSupportPoints][numComp]{};
+    real64 dGravHead_dP[numFluxSupportPoints]{};
+    real64 dGravHead_dC[numFluxSupportPoints][numComp]{};
+    PotGradZFormulation::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, useNewGravity, seri, sesri, sei, trans, dTrans_dPres, pres,
+                                                                   gravCoef, phaseVolFrac, dPhaseVolFrac, phaseMassDens, dPhaseMassDens,
+                                                                   phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, dPresGrad_dP,
+                                                                   dPresGrad_dC, dGravHead_dP, dGravHead_dC );
+
+    // *** upwinding ***
+
+    // choose upstream cell
+    k_up = (potGrad >= 0) ? 0 : 1;
+
+    localIndex const er_up  = seri[k_up];
+    localIndex const esr_up = sesri[k_up];
+    localIndex const ei_up  = sei[k_up];
+
+    real64 const mobility = phaseMob[er_up][esr_up][ei_up][ip];
+
+    // pressure gradient depends on all points in the stencil
+    for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+    {
+      dPhaseFlux_dP[ke] += dPresGrad_dP[ke] - dGravHead_dP[ke];
+      dPhaseFlux_dP[ke] *= mobility;
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        dPhaseFlux_dC[ke][jc] += dPresGrad_dC[ke][jc] - dGravHead_dC[ke][jc];
+        dPhaseFlux_dC[ke][jc] *= mobility;
+      }
+    }
+    // compute phase flux using upwind mobility.
+    phaseFlux = mobility * potGrad;
+
+    real64 const dMob_dP = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP];
+    arraySlice1d< real64 const, compflow::USD_PHASE_DC - 2 > dPhaseMobSub =
+      dPhaseMob[er_up][esr_up][ei_up][ip];
+
+    // add contribution from upstream cell mobility derivatives
+    dPhaseFlux_dP[k_up] += dMob_dP * potGrad;
+    for( integer jc = 0; jc < numComp; ++jc )
+    {
+      dPhaseFlux_dC[k_up][jc] += dPhaseMobSub[Deriv::dC+jc] * potGrad;
+    }
+
+    //distribute on phaseComponentFlux here
+    PhaseComponentFluxZFormulation::compute( ip, k_up, seri, sesri, sei, phaseCompFrac, dPhaseCompFrac, phaseFlux
+                                             , dPhaseFlux_dP, dPhaseFlux_dC, compFlux, dCompFlux_dP, dCompFlux_dC );
+
+  }
+};
+
+} // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities
+
+} // namespace geos
+
+
+#endif // GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PPUPHASEFLUXZFORMULATION_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseComponentFluxZFormulation.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseComponentFluxZFormulation.hpp
new file mode 100644
index 0000000000..16afaae3b3
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseComponentFluxZFormulation.hpp
@@ -0,0 +1,119 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file PhaseComponentFluxZFormulation.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PHASECOMPONENTFLUXZFORMULATION_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PHASECOMPONENTFLUXZFORMULATION_HPP
+
+#include "common/DataLayouts.hpp"
+#include "common/DataTypes.hpp"
+#include "constitutive/fluid/multifluid/Layouts.hpp"
+#include "mesh/ElementRegionManager.hpp"
+
+
+namespace geos
+{
+
+namespace isothermalCompositionalMultiphaseFVMKernelUtilities
+{
+
+template< typename VIEWTYPE >
+using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
+
+using Deriv = constitutive::multifluid::DerivativeOffset;
+
+struct PhaseComponentFluxZFormulation
+{
+  /**
+   * @brief Compute the component flux for a given phase
+   * @tparam numComp number of components
+   * @tparam numFluxSupportPoints number of flux support points
+   * @param ip phase index
+   * @param k_up uptream index for this phase
+   * @param seri arraySlice of the stencil-implied element region index
+   * @param sesri arraySlice of the stencil-implied element subregion index
+   * @param sei arraySlice of the stencil-implied element index
+   * @param phaseCompFrac phase component fraction
+   * @param dPhaseCompFrac derivative of phase component fraction wrt pressure, temperature, component fraction
+   * @param phaseFlux phase flux
+   * @param dPhaseFlux_dP derivative of phase flux wrt pressure
+   * @param dPhaseFlux_dC derivative of phase flux wrt comp density
+   * @param compFlux component flux
+   * @param dCompFlux_dP derivative of phase flux wrt pressure
+   * @param dCompFlux_dC derivative of phase flux wrt comp density
+   */
+  template< localIndex numComp, localIndex numFluxSupportPoints >
+  GEOS_HOST_DEVICE
+  static void
+  compute( localIndex const ip,
+           localIndex const k_up,
+           localIndex const ( &seri )[numFluxSupportPoints],
+           localIndex const ( &sesri )[numFluxSupportPoints],
+           localIndex const ( &sei )[numFluxSupportPoints],
+           ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac,
+           ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac,
+           real64 const & phaseFlux,
+           real64 const ( &dPhaseFlux_dP )[numFluxSupportPoints],
+           real64 const ( &dPhaseFlux_dC )[numFluxSupportPoints][numComp],
+           real64 ( & compFlux )[numComp],
+           real64 ( & dCompFlux_dP )[numFluxSupportPoints][numComp],
+           real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp] )
+  {
+    localIndex const er_up = seri[k_up];
+    localIndex const esr_up = sesri[k_up];
+    localIndex const ei_up = sei[k_up];
+
+    // slice some constitutive arrays to avoid too much indexing in component loop
+    arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE_COMP-3 > phaseCompFracSub =
+      phaseCompFrac[er_up][esr_up][ei_up][0][ip];
+    arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC-3 > dPhaseCompFracSub =
+      dPhaseCompFrac[er_up][esr_up][ei_up][0][ip];
+
+    // compute component fluxes and derivatives using upstream cell composition
+    for( integer ic = 0; ic < numComp; ++ic )
+    {
+      real64 const ycp = phaseCompFracSub[ic];
+      compFlux[ic] += phaseFlux * ycp;
+
+      // derivatives stemming from phase flux
+      for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+      {
+        dCompFlux_dP[ke][ic] += dPhaseFlux_dP[ke] * ycp;
+        for( integer jc = 0; jc < numComp; ++jc )
+        {
+          dCompFlux_dC[ke][ic][jc] += dPhaseFlux_dC[ke][jc] * ycp;
+        }
+      }
+
+      // additional derivatives stemming from upstream cell phase composition
+      dCompFlux_dP[k_up][ic] += phaseFlux * dPhaseCompFracSub[ic][Deriv::dP];
+
+      // convert derivatives of comp fraction w.r.t. comp fractions to derivatives w.r.t. comp densities
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        dCompFlux_dC[k_up][ic][jc] += phaseFlux * dPhaseCompFracSub[ic][Deriv::dC+jc];
+      }
+    }
+  }
+};
+
+} // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities
+
+} // namespace geos
+
+#endif // GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PHASECOMPONENTFLUXZFORMULATION_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseMobilityZFormulationKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseMobilityZFormulationKernel.hpp
new file mode 100644
index 0000000000..63f6f25d94
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseMobilityZFormulationKernel.hpp
@@ -0,0 +1,240 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file PhaseMobilityZFormulationKernel.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PHASEMOBILITYZFORMULATIONKERNEL_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PHASEMOBILITYZFORMULATIONKERNEL_HPP
+
+#include "physicsSolvers/fluidFlow/kernels/compositional/PropertyKernelBase.hpp"
+
+namespace geos
+{
+
+namespace isothermalCompositionalMultiphaseFVMKernels
+{
+
+/******************************** PhaseMobilityZFormulationKernel ********************************/
+
+/**
+ * @class PhaseMobilityZFormulationKernel
+ * @tparam NUM_COMP number of fluid components
+ * @tparam NUM_PHASE number of fluid phases
+ * @brief Defines the interface for the property kernel in charge of computing the phase mobilities
+ */
+template< integer NUM_COMP, integer NUM_PHASE >
+class PhaseMobilityZFormulationKernel : public isothermalCompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >
+{
+public:
+
+  using Base = isothermalCompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >;
+  using Base::numComp;
+
+  /// Compile time value for the number of phases
+  static constexpr integer numPhase = NUM_PHASE;
+
+  /**
+   * @brief Constructor
+   * @param[in] subRegion the element subregion
+   * @param[in] fluid the fluid model
+   * @param[in] relperm the relperm model
+   */
+  PhaseMobilityZFormulationKernel( ObjectManagerBase & subRegion,
+                                   constitutive::MultiFluidBase const & fluid,
+                                   constitutive::RelativePermeabilityBase const & relperm )
+    : Base(),
+    m_phaseVolFrac( subRegion.getField< fields::flow::phaseVolumeFraction >() ),
+    m_dPhaseVolFrac( subRegion.getField< fields::flow::dPhaseVolumeFraction >() ),
+    m_phaseDens( fluid.phaseDensity() ),
+    m_dPhaseDens( fluid.dPhaseDensity() ),
+    m_phaseVisc( fluid.phaseViscosity() ),
+    m_dPhaseVisc( fluid.dPhaseViscosity() ),
+    m_phaseRelPerm( relperm.phaseRelPerm() ),
+    m_dPhaseRelPerm_dPhaseVolFrac( relperm.dPhaseRelPerm_dPhaseVolFraction() ),
+    m_phaseMob( subRegion.getField< fields::flow::phaseMobility >() ),
+    m_dPhaseMob( subRegion.getField< fields::flow::dPhaseMobility >() )
+  {}
+
+  /**
+   * @brief Compute the phase mobilities in an element
+   * @tparam FUNC the type of the function that can be used to customize the kernel
+   * @param[in] ei the element index
+   * @param[in] phaseMobilityKernelOp the function used to customize the kernel
+   */
+  template< typename FUNC = NoOpFunc >
+  GEOS_HOST_DEVICE
+  void compute( localIndex const ei,
+                FUNC && phaseMobilityKernelOp = NoOpFunc{} ) const
+  {
+    using Deriv = constitutive::multifluid::DerivativeOffset;
+
+    arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const phaseDens = m_phaseDens[ei][0];
+    arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > const dPhaseDens = m_dPhaseDens[ei][0];
+    arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const phaseVisc = m_phaseVisc[ei][0];
+    arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > const dPhaseVisc = m_dPhaseVisc[ei][0];
+    arraySlice1d< real64 const, constitutive::relperm::USD_RELPERM - 2 > const phaseRelPerm = m_phaseRelPerm[ei][0];
+    arraySlice2d< real64 const, constitutive::relperm::USD_RELPERM_DS - 2 > const dPhaseRelPerm_dPhaseVolFrac = m_dPhaseRelPerm_dPhaseVolFrac[ei][0];
+    arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const phaseVolFrac = m_phaseVolFrac[ei];
+    arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const dPhaseVolFrac = m_dPhaseVolFrac[ei];
+    arraySlice1d< real64, compflow::USD_PHASE - 1 > const phaseMob = m_phaseMob[ei];
+    arraySlice2d< real64, compflow::USD_PHASE_DC - 1 > const dPhaseMob = m_dPhaseMob[ei];
+
+    real64 dRelPerm_dC[numComp]{};
+    real64 dDens_dC[numComp]{};
+    real64 dVisc_dC[numComp]{};
+
+    for( integer ip = 0; ip < numPhase; ++ip )
+    {
+
+      // compute the phase mobility only if the phase is present
+      bool const phaseExists = (phaseVolFrac[ip] > 0);
+      if( !phaseExists )
+      {
+        phaseMob[ip] = 0.0;
+        for( integer jc = 0; jc < numComp + 2; ++jc )
+        {
+          dPhaseMob[ip][jc] = 0.0;
+        }
+        continue;
+      }
+
+      real64 const density = phaseDens[ip];
+      real64 const dDens_dP = dPhaseDens[ip][Deriv::dP];
+      for( integer jc = 0; jc < numComp; ++jc )
+        dDens_dC[jc] = dPhaseDens[ip][Deriv::dC+jc];
+
+      real64 const viscosity = phaseVisc[ip];
+      real64 const dVisc_dP = dPhaseVisc[ip][Deriv::dP];
+      for( integer jc = 0; jc < numComp; ++jc )
+        dVisc_dC[jc] = dPhaseVisc[ip][Deriv::dC+jc];
+
+      real64 const relPerm = phaseRelPerm[ip];
+      real64 dRelPerm_dP = 0.0;
+      for( integer ic = 0; ic < numComp; ++ic )
+      {
+        dRelPerm_dC[ic] = 0.0;
+      }
+
+      for( integer jp = 0; jp < numPhase; ++jp )
+      {
+        real64 const dRelPerm_dS = dPhaseRelPerm_dPhaseVolFrac[ip][jp];
+        dRelPerm_dP += dRelPerm_dS * dPhaseVolFrac[jp][Deriv::dP];
+
+        for( integer jc = 0; jc < numComp; ++jc )
+        {
+          dRelPerm_dC[jc] += dRelPerm_dS * dPhaseVolFrac[jp][Deriv::dC+jc];
+        }
+      }
+
+      real64 const mobility = relPerm * density / viscosity;
+
+      phaseMob[ip] = mobility;
+      dPhaseMob[ip][Deriv::dP] = dRelPerm_dP * density / viscosity
+                                 + mobility * (dDens_dP / density - dVisc_dP / viscosity);
+
+      // compositional derivatives
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        dPhaseMob[ip][Deriv::dC+jc] = dRelPerm_dC[jc] * density / viscosity
+                                      + mobility * (dDens_dC[jc] / density - dVisc_dC[jc] / viscosity);
+      }
+
+      // call the lambda in the phase loop to allow the reuse of the relperm, density, viscosity, and mobility
+      // possible use: assemble the derivatives wrt temperature
+      phaseMobilityKernelOp( ip, phaseMob[ip], dPhaseMob[ip] );
+    }
+  }
+
+protected:
+
+  // inputs
+
+  /// Views on the phase volume fractions
+  arrayView2d< real64 const, compflow::USD_PHASE > m_phaseVolFrac;
+  arrayView3d< real64 const, compflow::USD_PHASE_DC > m_dPhaseVolFrac;
+
+  /// Views on the phase densities
+  arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > m_phaseDens;
+  arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > m_dPhaseDens;
+
+  /// Views on the phase viscosities
+  arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > m_phaseVisc;
+  arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > m_dPhaseVisc;
+
+  /// Views on the phase relative permeabilities
+  arrayView3d< real64 const, constitutive::relperm::USD_RELPERM > m_phaseRelPerm;
+  arrayView4d< real64 const, constitutive::relperm::USD_RELPERM_DS > m_dPhaseRelPerm_dPhaseVolFrac;
+
+  // outputs
+
+  /// Views on the phase mobilities
+  arrayView2d< real64, compflow::USD_PHASE > m_phaseMob;
+  arrayView3d< real64, compflow::USD_PHASE_DC > m_dPhaseMob;
+
+};
+
+/**
+ * @class PhaseMobilityZFormulationKernelFactory
+ */
+class PhaseMobilityZFormulationKernelFactory
+{
+public:
+
+  /**
+   * @brief Create a new kernel and launch
+   * @tparam POLICY the policy used in the RAJA kernel
+   * @param[in] numComp the number of fluid components
+   * @param[in] numPhase the number of fluid phases
+   * @param[in] subRegion the element subregion
+   * @param[in] fluid the fluid model
+   * @param[in] relperm the relperm model
+   */
+  template< typename POLICY >
+  static void
+  createAndLaunch( integer const numComp,
+                   integer const numPhase,
+                   ObjectManagerBase & subRegion,
+                   constitutive::MultiFluidBase const & fluid,
+                   constitutive::RelativePermeabilityBase const & relperm )
+  {
+    if( numPhase == 2 )
+    {
+      isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC )
+      {
+        integer constexpr NUM_COMP = NC();
+        PhaseMobilityZFormulationKernel< NUM_COMP, 2 > kernel( subRegion, fluid, relperm );
+        PhaseMobilityZFormulationKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel );
+      } );
+    }
+    else if( numPhase == 3 )
+    {
+      isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC )
+      {
+        integer constexpr NUM_COMP = NC();
+        PhaseMobilityZFormulationKernel< NUM_COMP, 3 > kernel( subRegion, fluid, relperm );
+        PhaseMobilityZFormulationKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel );
+      } );
+    }
+  }
+};
+
+} // namespace isothermalCompositionalMultiphaseFVMKernels
+
+} // namespace geos
+
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PHASEMOBILITYZFORMULATIONKERNEL_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseVolumeFractionZFormulationKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseVolumeFractionZFormulationKernel.hpp
new file mode 100644
index 0000000000..819d21c8ab
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PhaseVolumeFractionZFormulationKernel.hpp
@@ -0,0 +1,230 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file PhaseVolumeFractionZFormulationKernel.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PHASEVOLUMEFRACTIONZFORMULATIONKERNEL_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PHASEVOLUMEFRACTIONZFORMULATIONKERNEL_HPP
+
+#include "physicsSolvers/fluidFlow/kernels/compositional/PropertyKernelBase.hpp"
+
+namespace geos
+{
+
+namespace isothermalCompositionalMultiphaseBaseKernels
+{
+
+/******************************** PhaseVolumeFractionZFormulationKernel ********************************/
+
+/**
+ * @class PhaseVolumeFractionZFormulationKernel
+ * @tparam NUM_COMP number of fluid components
+ * @tparam NUM_PHASE number of fluid phases
+ * @brief Define the interface for the property kernel in charge of computing the phase volume fractions
+ */
+template< integer NUM_COMP, integer NUM_PHASE >
+class PhaseVolumeFractionZFormulationKernel : public PropertyKernelBase< NUM_COMP >
+{
+public:
+
+  using Base = PropertyKernelBase< NUM_COMP >;
+  using Base::numComp;
+
+  /// Compile time value for the number of phases
+  static constexpr integer numPhase = NUM_PHASE;
+
+  /**
+   * @brief Constructor
+   * @param[in] subRegion the element subregion
+   * @param[in] fluid the fluid model
+   */
+  PhaseVolumeFractionZFormulationKernel( ObjectManagerBase & subRegion,
+                                         constitutive::MultiFluidBase const & fluid )
+    : Base(),
+    m_phaseVolFrac( subRegion.getField< fields::flow::phaseVolumeFraction >() ),
+    m_dPhaseVolFrac( subRegion.getField< fields::flow::dPhaseVolumeFraction >() ),
+    m_phaseFrac( fluid.phaseFraction() ),
+    m_dPhaseFrac( fluid.dPhaseFraction() ),
+    m_phaseDens( fluid.phaseDensity() ),
+    m_dPhaseDens( fluid.dPhaseDensity() ),
+    m_totalDens( fluid.totalDensity() ),
+    m_dTotalDens( fluid.dTotalDensity() )
+  {}
+
+  /**
+   * @brief Compute the phase volume fractions in an element
+   * @tparam FUNC the type of the function that can be used to customize the kernel
+   * @param[in] ei the element index
+   * @param[in] phaseVolFractionKernelOp the function used to customize the kernel
+   */
+  template< typename FUNC = NoOpFunc >
+  GEOS_HOST_DEVICE
+  real64 compute( localIndex const ei ) const
+  {
+    using Deriv = constitutive::multifluid::DerivativeOffset;
+
+    arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const phaseFrac = m_phaseFrac[ei][0];
+    arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > const dPhaseFrac = m_dPhaseFrac[ei][0];
+    arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const phaseDens = m_phaseDens[ei][0];
+    arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > const dPhaseDens = m_dPhaseDens[ei][0];
+    real64 const totalDensity = m_totalDens[ei][0];
+    arraySlice1d< real64 const, constitutive::multifluid::USD_FLUID_DC - 2 > const dTotalDens = m_dTotalDens[ei][0];
+    arraySlice1d< real64, compflow::USD_PHASE - 1 > const phaseVolFrac = m_phaseVolFrac[ei];
+    arraySlice2d< real64, compflow::USD_PHASE_DC - 1 > const dPhaseVolFrac = m_dPhaseVolFrac[ei];
+
+    real64 maxDeltaPhaseVolFrac = 0.0;
+
+    for( integer ip = 0; ip < numPhase; ++ip )
+    {
+
+      // set the saturation to zero if the phase is absent
+      bool const phaseExists = (phaseFrac[ip] > 0);
+      if( !phaseExists )
+      {
+        phaseVolFrac[ip] = 0.;
+        for( integer jc = 0; jc < numComp+2; ++jc )
+        {
+          dPhaseVolFrac[ip][jc] = 0.;
+        }
+        continue;
+      }
+
+      // store old saturation to compute change later
+      real64 const satOld = phaseVolFrac[ip];
+
+      // compute saturation and derivatives  S_p = rho_t * (nu_p / rho_p)
+      phaseVolFrac[ip] = totalDensity * phaseFrac[ip] / phaseDens[ip];
+
+      dPhaseVolFrac[ip][Deriv::dP] = phaseVolFrac[ip] *
+                                     (dTotalDens[Deriv::dP] / totalDensity +  dPhaseFrac[ip][Deriv::dP] / phaseFrac[ip] - dPhaseDens[ip][Deriv::dP] / phaseDens[ip]);
+
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        dPhaseVolFrac[ip][Deriv::dC+jc] = phaseVolFrac[ip] *
+                                          (dTotalDens[Deriv::dC+jc] / totalDensity +  dPhaseFrac[ip][Deriv::dC+jc] / phaseFrac[ip] - dPhaseDens[ip][Deriv::dC+jc] / phaseDens[ip]);
+      }
+
+      // call the lambda in the phase loop to allow the reuse of the phaseVolFrac and totalDensity
+      // possible use: assemble the derivatives wrt temperature
+      //phaseVolFractionKernelOp( ip, phaseVolFrac[ip], phaseDensInv, totalDensity );
+
+      real64 const deltaPhaseVolFrac = LvArray::math::abs( phaseVolFrac[ip] - satOld );
+
+      if( maxDeltaPhaseVolFrac < deltaPhaseVolFrac )
+      {
+        maxDeltaPhaseVolFrac = deltaPhaseVolFrac;
+      }
+    }
+    return maxDeltaPhaseVolFrac;
+  }
+
+  /**
+   * @brief Performs the kernel launch
+   * @tparam POLICY the policy used in the RAJA kernels
+   * @tparam KERNEL_TYPE the kernel type
+   * @param[in] numElems the number of elements
+   * @param[inout] kernelComponent the kernel component providing access to the compute function
+   */
+  template< typename POLICY, typename KERNEL_TYPE >
+  static real64
+  launch( localIndex const numElems,
+          KERNEL_TYPE const & kernelComponent )
+  {
+    RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaPhaseVolFrac( 0.0 );
+    forAll< POLICY >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei )
+    {
+      real64 const deltaPhaseVolFrac = kernelComponent.compute( ei );
+      maxDeltaPhaseVolFrac.max( deltaPhaseVolFrac );
+    } );
+    return maxDeltaPhaseVolFrac.get();
+  }
+
+protected:
+
+  // outputs
+
+  /// Views on phase volume fractions
+  arrayView2d< real64, compflow::USD_PHASE > m_phaseVolFrac;
+  arrayView3d< real64, compflow::USD_PHASE_DC > m_dPhaseVolFrac;
+
+  // inputs
+
+  /// Views on phase fractions
+  arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > m_phaseFrac;
+  arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > m_dPhaseFrac;
+
+  /// Views on phase densities
+  arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > m_phaseDens;
+  arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > m_dPhaseDens;
+
+  /// Views on the total density
+  arrayView2d< real64 const, constitutive::multifluid::USD_FLUID > const m_totalDens;
+  arrayView3d< real64 const, constitutive::multifluid::USD_FLUID_DC > const m_dTotalDens;
+
+};
+
+/**
+ * @class PhaseVolumeFractionZFormulationKernelFactory
+ */
+class PhaseVolumeFractionZFormulationKernelFactory
+{
+public:
+
+  /**
+   * @brief Create a new kernel and launch
+   * @tparam POLICY the policy used in the RAJA kernel
+   * @param[in] numComp the number of fluid components
+   * @param[in] numPhase the number of fluid phases
+   * @param[in] subRegion the element subregion
+   * @param[in] fluid the fluid model
+   */
+  template< typename POLICY >
+  static real64
+  createAndLaunch( integer const numComp,
+                   integer const numPhase,
+                   ObjectManagerBase & subRegion,
+                   constitutive::MultiFluidBase const & fluid )
+  {
+    real64 maxDeltaPhaseVolFrac = 0.0;
+    if( numPhase == 2 )
+    {
+      internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC )
+      {
+        integer constexpr NUM_COMP = NC();
+        PhaseVolumeFractionZFormulationKernel< NUM_COMP, 2 > kernel( subRegion, fluid );
+        maxDeltaPhaseVolFrac = PhaseVolumeFractionZFormulationKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel );
+      } );
+    }
+    else if( numPhase == 3 )
+    {
+      internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC )
+      {
+        integer constexpr NUM_COMP = NC();
+        PhaseVolumeFractionZFormulationKernel< NUM_COMP, 3 > kernel( subRegion, fluid );
+        maxDeltaPhaseVolFrac = PhaseVolumeFractionZFormulationKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel );
+      } );
+    }
+    return maxDeltaPhaseVolFrac;
+  }
+};
+
+} // namespace isothermalCompositionalMultiphaseBaseKernels
+
+} // namespace geos
+
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_PHASEVOLUMEFRACTIONZFORMULATIONKERNEL_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PotGradZFormulation.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PotGradZFormulation.hpp
new file mode 100644
index 0000000000..4cd4a758ea
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/PotGradZFormulation.hpp
@@ -0,0 +1,216 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file PotGradZFormulation.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_POTGRADZFORMULATION_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_POTGRADZFORMULATION_HPP
+
+#include "common/DataLayouts.hpp"
+#include "common/DataTypes.hpp"
+#include "constitutive/fluid/multifluid/Layouts.hpp"
+#include "constitutive/capillaryPressure/layouts.hpp"
+#include "mesh/ElementRegionManager.hpp"
+
+
+namespace geos
+{
+
+namespace isothermalCompositionalMultiphaseFVMKernelUtilities
+{
+
+template< typename VIEWTYPE >
+using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
+
+using Deriv = constitutive::multifluid::DerivativeOffset;
+
+struct PotGradZFormulation
+{
+  template< integer numComp, integer numFluxSupportPoints >
+  GEOS_HOST_DEVICE
+  static void
+  compute ( integer const numPhase,
+            integer const ip,
+            integer const hasCapPressure,
+            integer const useNewGravity,
+            localIndex const ( &seri )[numFluxSupportPoints],
+            localIndex const ( &sesri )[numFluxSupportPoints],
+            localIndex const ( &sei )[numFluxSupportPoints],
+            real64 const ( &trans )[numFluxSupportPoints],
+            real64 const ( &dTrans_dPres )[numFluxSupportPoints],
+            ElementViewConst< arrayView1d< real64 const > > const & pres,
+            ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
+            ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
+            ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
+            ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
+            ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
+            ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
+            ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
+            real64 & potGrad,
+            real64 ( & dPresGrad_dP )[numFluxSupportPoints],
+            real64 ( & dPresGrad_dC )[numFluxSupportPoints][numComp],
+            real64 ( & dGravHead_dP )[numFluxSupportPoints],
+            real64 ( & dGravHead_dC )[numFluxSupportPoints][numComp] )
+  {
+    // assign derivatives arrays to zero
+    for( integer i = 0; i < numFluxSupportPoints; ++i )
+    {
+      dPresGrad_dP[i] = 0.0;
+      dGravHead_dP[i] = 0.0;
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        dPresGrad_dC[i][jc] = 0.0;
+        dGravHead_dC[i][jc] = 0.0;
+      }
+    }
+
+    // create local work arrays
+    real64 densMean = 0.0;
+    real64 dDensMean_dP[numFluxSupportPoints]{};
+    real64 dDensMean_dC[numFluxSupportPoints][numComp]{};
+
+    real64 presGrad = 0.0;
+    real64 gravHead = 0.0;
+    real64 dCapPressure_dC[numComp]{};
+
+    calculateMeanDensity( useNewGravity, ip, seri, sesri, sei, phaseVolFrac, phaseMassDens, dPhaseMassDens, densMean, dDensMean_dP, dDensMean_dC );
+
+    /// compute the TPFA potential difference
+    for( integer i = 0; i < numFluxSupportPoints; i++ )
+    {
+      localIndex const er  = seri[i];
+      localIndex const esr = sesri[i];
+      localIndex const ei  = sei[i];
+
+      // capillary pressure
+      real64 capPressure     = 0.0;
+      real64 dCapPressure_dP = 0.0;
+
+      for( integer ic = 0; ic < numComp; ++ic )
+      {
+        dCapPressure_dC[ic] = 0.0;
+      }
+
+      if( hasCapPressure )
+      {
+        capPressure = phaseCapPressure[er][esr][ei][0][ip];
+
+        for( integer jp = 0; jp < numPhase; ++jp )
+        {
+          real64 const dCapPressure_dS = dPhaseCapPressure_dPhaseVolFrac[er][esr][ei][0][ip][jp];
+          dCapPressure_dP += dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dP];
+
+          for( integer jc = 0; jc < numComp; ++jc )
+          {
+            dCapPressure_dC[jc] += dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dC+jc];
+          }
+        }
+      }
+
+      presGrad += trans[i] * (pres[er][esr][ei] - capPressure);
+      dPresGrad_dP[i] += trans[i] * (1 - dCapPressure_dP)
+                         + dTrans_dPres[i] * (pres[er][esr][ei] - capPressure);
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        dPresGrad_dC[i][jc] += -trans[i] * dCapPressure_dC[jc];
+      }
+
+      real64 const gravD     = trans[i] * gravCoef[er][esr][ei];
+      real64 const dGravD_dP = dTrans_dPres[i] * gravCoef[er][esr][ei];
+
+      // the density used in the potential difference is always a mass density
+      // unlike the density used in the phase mobility, which is a mass density
+      // if useMass == 1 and a molar density otherwise
+      gravHead += densMean * gravD;
+
+      // need to add contributions from both cells the mean density depends on
+      for( integer j = 0; j < numFluxSupportPoints; ++j )
+      {
+        dGravHead_dP[j] += dDensMean_dP[j] * gravD + dGravD_dP * densMean;
+        for( integer jc = 0; jc < numComp; ++jc )
+        {
+          dGravHead_dC[j][jc] += dDensMean_dC[j][jc] * gravD;
+        }
+      }
+    }
+
+    // compute phase potential gradient
+    potGrad = presGrad - gravHead;
+
+  }
+
+  template< integer numComp, integer numFluxSupportPoints >
+  GEOS_HOST_DEVICE
+  static void
+  calculateMeanDensity( integer const useNewGravity,
+                        integer const ip,
+                        localIndex const ( &seri )[numFluxSupportPoints],
+                        localIndex const ( &sesri )[numFluxSupportPoints],
+                        localIndex const ( &sei )[numFluxSupportPoints],
+                        ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
+                        ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
+                        ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
+                        real64 & densMean, real64 ( & dDensMean_dP)[numFluxSupportPoints], real64 ( & dDensMean_dC )[numFluxSupportPoints][numComp] )
+  {
+    // calculate quantities on primary connected cells
+    integer denom = 0;
+    for( integer i = 0; i < numFluxSupportPoints; ++i )
+    {
+      localIndex const er  = seri[i];
+      localIndex const esr = sesri[i];
+      localIndex const ei  = sei[i];
+
+      bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0);
+      if( useNewGravity && !phaseExists )
+      {
+        continue;
+      }
+
+      // density
+      real64 const density  = phaseMassDens[er][esr][ei][0][ip];
+      real64 const dDens_dP = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP];
+
+      // average density and derivatives
+      densMean += density;
+      dDensMean_dP[i] = dDens_dP;
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        dDensMean_dC[i][jc] = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dC+jc];
+      }
+      denom++;
+    }
+    if( denom > 1 )
+    {
+      densMean /= denom;
+      for( integer i = 0; i < numFluxSupportPoints; ++i )
+      {
+        dDensMean_dP[i] /= denom;
+        for( integer jc = 0; jc < numComp; ++jc )
+        {
+          dDensMean_dC[i][jc] /= denom;
+        }
+      }
+    }
+  }
+
+};
+
+} // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities
+
+} // namespace geos
+
+#endif // GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_POTGRADZFORMULATION_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingAndCheckingZFormulationKernelBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingAndCheckingZFormulationKernelBase.hpp
new file mode 100644
index 0000000000..3b3d7f43d4
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingAndCheckingZFormulationKernelBase.hpp
@@ -0,0 +1,182 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file SolutionScalingAndCheckingZFormulationKernelBase.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_SOLUTIONSCALINGANDCHECKINGZFORMULATIONKERNELBASE_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_SOLUTIONSCALINGANDCHECKINGZFORMULATIONKERNELBASE_HPP
+
+#include "common/DataLayouts.hpp"
+#include "common/DataTypes.hpp"
+#include "common/GEOS_RAJA_Interface.hpp"
+#include "mesh/ElementSubRegionBase.hpp"
+
+namespace geos
+{
+
+namespace isothermalCompositionalMultiphaseBaseKernels
+{
+
+/**
+ * @class SolutionScalingAndCheckingZFormulationKernelBase
+ * @brief Define the kernel for scaling the solution and check its validity
+ */
+template< typename TYPE >
+class SolutionScalingAndCheckingZFormulationKernelBase
+{
+public:
+
+  /**
+   * @brief Create a new kernel instance
+   * @param[in] rankOffset the rank offset
+   * @param[in] numComp the number of components
+   * @param[in] dofKey the dof key to get dof numbers
+   * @param[in] subRegion the subRegion
+   * @param[in] localSolution the Newton update
+   * @param[in] pressure the pressure vector
+   * @param[in] compFrac the component density vector
+   * @param[in] pressureScalingFactor the pressure local scaling factor
+   * @param[in] compFracScalingFactor the component local scaling factor
+   */
+  SolutionScalingAndCheckingZFormulationKernelBase( globalIndex const rankOffset,
+                                                    integer const numComp,
+                                                    string const dofKey,
+                                                    ElementSubRegionBase const & subRegion,
+                                                    arrayView1d< real64 const > const localSolution,
+                                                    arrayView1d< real64 const > const pressure,
+                                                    arrayView2d< real64 const, compflow::USD_COMP > const compFrac,
+                                                    arrayView1d< real64 > pressureScalingFactor,
+                                                    arrayView1d< real64 > compFracScalingFactor )
+    : m_rankOffset( rankOffset ),
+    m_numComp( numComp ),
+    m_dofNumber( subRegion.getReference< array1d< globalIndex > >( dofKey ) ),
+    m_ghostRank( subRegion.ghostRank() ),
+    m_localSolution( localSolution ),
+    m_pressure( pressure ),   // not passed with fields::flow to be able to reuse this for wells
+    m_compFrac( compFrac ),   // same here
+    m_pressureScalingFactor( pressureScalingFactor ),
+    m_compFracScalingFactor( compFracScalingFactor )
+  { }
+
+  /**
+   * @struct StackVariables
+   * @brief Kernel variables located on the stack
+   */
+  struct StackVariables
+  {
+    GEOS_HOST_DEVICE
+    StackVariables()
+    { }
+
+    StackVariables( real64 _localMinVal )
+      :
+      localMinVal( _localMinVal )
+    { }
+
+    /// Index of the local row corresponding to this element
+    localIndex localRow;
+
+    /// The local value
+    TYPE localMinVal;
+  };
+
+  /**
+   * @brief Performs the setup phase for the kernel.
+   * @param[in] ei the element index
+   * @param[in] stack the stack variables
+   */
+  GEOS_HOST_DEVICE
+  void setup( localIndex const ei,
+              StackVariables & stack ) const
+  {
+    stack.localMinVal = 1;
+
+    // set row index and degrees of freedom indices for this element
+    stack.localRow = m_dofNumber[ei] - m_rankOffset;
+  }
+
+  /**
+   * @brief Getter for the ghost rank
+   * @param[in] i the looping index of the element/node/face
+   * @return the ghost rank of the element/node/face
+   */
+  GEOS_HOST_DEVICE
+  integer ghostRank( localIndex const i ) const
+  { return m_ghostRank( i ); }
+
+  /**
+   * @brief Performs the kernel launch
+   * @tparam POLICY the policy used in the RAJA kernels
+   * @tparam KERNEL_TYPE the kernel type
+   * @param[in] numElems the number of elements
+   * @param[inout] kernelComponent the kernel component providing access to the compute function
+   */
+  template< typename POLICY, typename KERNEL_TYPE >
+  static TYPE
+  launch( localIndex const numElems,
+          KERNEL_TYPE const & kernelComponent )
+  {
+    RAJA::ReduceMin< ReducePolicy< POLICY >, TYPE > minVal( 1 );
+    forAll< POLICY >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei )
+    {
+      if( kernelComponent.ghostRank( ei ) >= 0 )
+      {
+        return;
+      }
+
+      StackVariables stack;
+      kernelComponent.setup( ei, stack );
+      kernelComponent.compute( ei, stack );
+      minVal.min( stack.localMinVal );
+    } );
+
+    return minVal.get();
+  }
+
+protected:
+
+  /// Offset for my MPI rank
+  globalIndex const m_rankOffset;
+
+  /// Number of components
+  real64 const m_numComp;
+
+  /// View on the dof numbers
+  arrayView1d< globalIndex const > const m_dofNumber;
+
+  /// View on the ghost ranks
+  arrayView1d< integer const > const m_ghostRank;
+
+  /// View on the local residual
+  arrayView1d< real64 const > const m_localSolution;
+
+  /// View on the primary variables
+  arrayView1d< real64 const > const m_pressure;
+  arrayView2d< real64 const, compflow::USD_COMP > const m_compFrac;
+
+  /// View on the scaling factors
+  arrayView1d< real64 > const m_pressureScalingFactor;
+  arrayView1d< real64 > const m_compFracScalingFactor;
+
+};
+
+} // namespace isothermalCompositionalMultiphaseBaseKernels
+
+} // namespace geos
+
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_SOLUTIONSCALINGANDCHECKINGZFORMULATIONKERNELBASE_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingZFormulationKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingZFormulationKernel.hpp
new file mode 100644
index 0000000000..4786871950
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingZFormulationKernel.hpp
@@ -0,0 +1,348 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file SolutionScalingZFormulationKernel.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_SOLUTIONSCALINGZFORMULATIONKERNEL_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_SOLUTIONSCALINGZFORMULATIONKERNEL_HPP
+
+#include "physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingAndCheckingZFormulationKernelBase.hpp"
+#include "physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp"
+
+namespace geos
+{
+
+namespace isothermalCompositionalMultiphaseBaseKernels
+{
+
+/******************************** SolutionScalingZFormulationKernel ********************************/
+
+/**
+ * @class SolutionScalingZFormulationKernel
+ * @brief Define the kernel for scaling the Newton update
+ */
+class SolutionScalingZFormulationKernel : public SolutionScalingAndCheckingZFormulationKernelBase< real64 >
+{
+public:
+
+  using Base = SolutionScalingAndCheckingZFormulationKernelBase< real64 >;
+  using Base::m_rankOffset;
+  using Base::m_numComp;
+  using Base::m_dofNumber;
+  using Base::m_ghostRank;
+  using Base::m_localSolution;
+  using Base::m_pressure;
+  using Base::m_compFrac;
+  using Base::m_pressureScalingFactor;
+  using Base::m_compFracScalingFactor;
+
+  /**
+   * @brief Create a new kernel instance
+   * @param[in] maxRelativePresChange the max allowed relative pressure change
+   * @param[in] maxAbsolutePresChange the max allowed absolute pressure change
+   * @param[in] maxCompFracChange the max allowed comp fraction change
+   * @param[in] rankOffset the rank offset
+   * @param[in] numComp the number of components
+   * @param[in] dofKey the dof key to get dof numbers
+   * @param[in] subRegion the subRegion
+   * @param[in] localSolution the Newton update
+   * @param[in] pressure the pressure vector
+   * @param[in] compFrac the component density vector
+   * @param[in] pressureScalingFactor the pressure local scaling factor
+   * @param[in] compFracScalingFactor the component density local scaling factor
+   */
+  SolutionScalingZFormulationKernel( real64 const maxRelativePresChange,
+                                     real64 const maxAbsolutePresChange,
+                                     real64 const maxCompFracChange,
+                                     globalIndex const rankOffset,
+                                     integer const numComp,
+                                     string const dofKey,
+                                     ElementSubRegionBase const & subRegion,
+                                     arrayView1d< real64 const > const localSolution,
+                                     arrayView1d< real64 const > const pressure,
+                                     arrayView2d< real64 const, compflow::USD_COMP > const compFrac,
+                                     arrayView1d< real64 > pressureScalingFactor,
+                                     arrayView1d< real64 > compFracScalingFactor )
+    : Base( rankOffset,
+            numComp,
+            dofKey,
+            subRegion,
+            localSolution,
+            pressure,
+            compFrac,
+            pressureScalingFactor,
+            compFracScalingFactor ),
+    m_maxRelativePresChange( maxRelativePresChange ),
+    m_maxAbsolutePresChange( maxAbsolutePresChange ),
+    m_maxCompFracChange( maxCompFracChange )
+  {}
+
+  /**
+   * @struct StackVariables
+   * @brief Kernel variables located on the stack
+   */
+  struct StackVariables : public Base::StackVariables
+  {
+    GEOS_HOST_DEVICE
+    StackVariables()
+    { }
+
+    StackVariables( real64 _localMinVal,
+                    real64 _localMaxDeltaPres,
+                    real64 _localMaxDeltaTemp,
+                    real64 _localMaxDeltaCompFrac,
+                    real64 _localMinPresScalingFactor,
+                    real64 _localMinTempScalingFactor,
+                    real64 _localMinCompFracScalingFactor )
+      :
+      Base::StackVariables( _localMinVal ),
+      localMaxDeltaPres( _localMaxDeltaPres ),
+      localMaxDeltaTemp( _localMaxDeltaTemp ),
+      localMaxDeltaCompFrac( _localMaxDeltaCompFrac ),
+      localMinPresScalingFactor( _localMinPresScalingFactor ),
+      localMinTempScalingFactor( _localMinTempScalingFactor ),
+      localMinCompFracScalingFactor( _localMinCompFracScalingFactor )
+    { }
+
+    real64 localMaxDeltaPres;
+    real64 localMaxDeltaTemp;
+    real64 localMaxDeltaCompFrac;
+
+    real64 localMinPresScalingFactor;
+    real64 localMinTempScalingFactor;
+    real64 localMinCompFracScalingFactor;
+
+  };
+
+  /**
+   * @brief Performs the kernel launch
+   * @tparam POLICY the policy used in the RAJA kernels
+   * @tparam KERNEL_TYPE the kernel type
+   * @param[in] numElems the number of elements
+   * @param[inout] kernelComponent the kernel component providing access to the compute function
+   */
+  template< typename POLICY, typename KERNEL_TYPE >
+  static StackVariables
+  launch( localIndex const numElems,
+          KERNEL_TYPE const & kernelComponent )
+  {
+    RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > globalScalingFactor( 1.0 );
+
+    RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaPres( 0.0 );
+    RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaTemp( 0.0 );
+    RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaCompFrac( 0.0 );
+
+    RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > minPresScalingFactor( 1.0 );
+    RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > minTempScalingFactor( 1.0 );
+    RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > minCompFracScalingFactor( 1.0 );
+
+    forAll< POLICY >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei )
+    {
+      if( kernelComponent.ghostRank( ei ) >= 0 )
+      {
+        return;
+      }
+
+      StackVariables stack;
+      kernelComponent.setup( ei, stack );
+      kernelComponent.compute( ei, stack );
+
+      globalScalingFactor.min( stack.localMinVal );
+
+      maxDeltaPres.max( stack.localMaxDeltaPres );
+      maxDeltaTemp.max( stack.localMaxDeltaTemp );
+      maxDeltaCompFrac.max( stack.localMaxDeltaCompFrac );
+
+      minPresScalingFactor.min( stack.localMinPresScalingFactor );
+      minTempScalingFactor.min( stack.localMinTempScalingFactor );
+      minCompFracScalingFactor.min( stack.localMinCompFracScalingFactor );
+    } );
+
+    return StackVariables( globalScalingFactor.get(),
+                           maxDeltaPres.get(),
+                           maxDeltaTemp.get(),
+                           maxDeltaCompFrac.get(),
+                           minPresScalingFactor.get(),
+                           minTempScalingFactor.get(),
+                           minCompFracScalingFactor.get() );
+  }
+
+  GEOS_HOST_DEVICE
+  void setup( localIndex const ei,
+              StackVariables & stack ) const
+  {
+    Base::setup( ei, stack );
+
+    stack.localMaxDeltaPres = 0.0;
+    stack.localMaxDeltaTemp = 0.0;
+    stack.localMaxDeltaCompFrac = 0.0;
+
+    stack.localMinPresScalingFactor = 1.0;
+    stack.localMinTempScalingFactor = 1.0;
+    stack.localMinCompFracScalingFactor = 1.0;
+  }
+
+  /**
+   * @brief Compute the local value
+   * @param[in] ei the element index
+   * @param[inout] stack the stack variables
+   */
+  GEOS_HOST_DEVICE
+  void compute( localIndex const ei,
+                StackVariables & stack ) const
+  {
+    computeScalingFactor( ei, stack );
+  }
+
+  /**
+   * @brief Compute the local value of the scaling factor
+   * @tparam FUNC the type of the function that can be used to customize the kernel
+   * @param[in] ei the element index
+   * @param[inout] stack the stack variables
+   * @param[in] kernelOp the function used to customize the kernel
+   */
+  template< typename FUNC = NoOpFunc >
+  GEOS_HOST_DEVICE
+  void computeScalingFactor( localIndex const ei,
+                             StackVariables & stack,
+                             FUNC && kernelOp = NoOpFunc{} ) const
+  {
+    real64 constexpr eps = minDensForDivision;
+
+    // compute the change in pressure
+    real64 const pres = m_pressure[ei];
+    real64 const absPresChange = LvArray::math::abs( m_localSolution[stack.localRow] );
+    if( stack.localMaxDeltaPres < absPresChange )
+    {
+      stack.localMaxDeltaPres = absPresChange;
+    }
+
+    // compute pressure scaling factor
+    real64 presScalingFactor = 1.0;
+    // when enabled, absolute change scaling has a priority over relative change
+    if( m_maxAbsolutePresChange > 0.0 ) // maxAbsolutePresChange <= 0.0 means that absolute scaling is disabled
+    {
+      if( absPresChange > m_maxAbsolutePresChange )
+      {
+        presScalingFactor = m_maxAbsolutePresChange / absPresChange;
+      }
+    }
+    else if( pres > eps )
+    {
+      real64 const relativePresChange = absPresChange / pres;
+      if( relativePresChange > m_maxRelativePresChange )
+      {
+        presScalingFactor = m_maxRelativePresChange / relativePresChange;
+      }
+    }
+    m_pressureScalingFactor[ei] = presScalingFactor;
+    if( stack.localMinVal > presScalingFactor )
+    {
+      stack.localMinVal = presScalingFactor;
+    }
+    if( stack.localMinPresScalingFactor > presScalingFactor )
+    {
+      stack.localMinPresScalingFactor = presScalingFactor;
+    }
+
+    // Component Fraction Change Scaling
+    m_compFracScalingFactor[ei] = 1.0;
+
+    for( integer ic = 0; ic < m_numComp; ++ic )
+    {
+      // compute scaling factor based on relative change in component densities
+      real64 const absCompFracChange = LvArray::math::abs( m_localSolution[stack.localRow + ic + 1] );
+      if( stack.localMaxDeltaCompFrac < absCompFracChange )
+      {
+        stack.localMaxDeltaCompFrac = absCompFracChange;
+      }
+
+      if( absCompFracChange > m_maxCompFracChange && absCompFracChange > eps )
+      {
+        real64 const compScalingFactor = m_maxCompFracChange / absCompFracChange;
+        m_compFracScalingFactor[ei] = LvArray::math::min( m_compFracScalingFactor[ei], compScalingFactor );
+        if( stack.localMinVal > compScalingFactor )
+        {
+          stack.localMinVal = compScalingFactor;
+        }
+        if( stack.localMinCompFracScalingFactor > compScalingFactor )
+        {
+          stack.localMinCompFracScalingFactor = compScalingFactor;
+        }
+      }
+    }
+
+    // compute the scaling factor for other vars, such as temperature
+    kernelOp();
+  }
+
+protected:
+
+  /// Max allowed changes in primary variables
+  real64 const m_maxRelativePresChange;
+  real64 const m_maxAbsolutePresChange;
+  real64 const m_maxCompFracChange;
+
+};
+
+/**
+ * @class SolutionScalingZFormulationKernelFactory
+ */
+class SolutionScalingZFormulationKernelFactory
+{
+public:
+
+  /*
+   * @brief Create and launch the kernel computing the scaling factor
+   * @tparam POLICY the kernel policy
+   * @param[in] maxRelativePresChange the max allowed relative pressure change
+   * @param[in] maxAbsolutePresChange the max allowed absolute pressure change
+   * @param[in] maxCompFracChange the max allowed comp fraction change
+   * @param[in] rankOffset the rank offset
+   * @param[in] numComp the number of components
+   * @param[in] dofKey the dof key to get dof numbers
+   * @param[in] subRegion the subRegion
+   * @param[in] localSolution the Newton update
+   * @return the scaling factor
+   */
+  template< typename POLICY >
+  static SolutionScalingZFormulationKernel::StackVariables
+  createAndLaunch( real64 const maxRelativePresChange,
+                   real64 const maxAbsolutePresChange,
+                   real64 const maxCompFracChange,
+                   arrayView1d< real64 const > const pressure,
+                   arrayView2d< real64 const, compflow::USD_COMP > const compFrac,
+                   arrayView1d< real64 > pressureScalingFactor,
+                   arrayView1d< real64 > compFracScalingFactor,
+                   globalIndex const rankOffset,
+                   integer const numComp,
+                   string const dofKey,
+                   ElementSubRegionBase & subRegion,
+                   arrayView1d< real64 const > const localSolution )
+  {
+    SolutionScalingZFormulationKernel kernel( maxRelativePresChange, maxAbsolutePresChange, maxCompFracChange, rankOffset,
+                                              numComp, dofKey, subRegion, localSolution, pressure, compFrac, pressureScalingFactor, compFracScalingFactor );
+    return SolutionScalingZFormulationKernel::launch< POLICY >( subRegion.size(), kernel );
+  }
+};
+
+} // namespace isothermalCompositionalMultiphaseBaseKernels
+
+} // namespace geos
+
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_SOLUTIONSCALINGZFORMULATIONKERNEL_HPP