diff --git a/modules/navier_stokes/doc/content/bib/navier_stokes.bib b/modules/navier_stokes/doc/content/bib/navier_stokes.bib
index 1990737d60b8..4740ab463f2a 100644
--- a/modules/navier_stokes/doc/content/bib/navier_stokes.bib
+++ b/modules/navier_stokes/doc/content/bib/navier_stokes.bib
@@ -336,3 +336,44 @@ @phdthesis{jasak1996error
   school={Imperial College London (University of London)}
 }
 
+@InProceedings{suikkanen,
+      author = {H. Suikkanen and V. Rintala and R. Kyrki-Rajamaki},
+       title = {{Development of a Coupled Multi-Physics Code System for Pebble Bed Reactor Core Modeling}},
+  booktitle  = {{Proceedings of the HTR 2014}},
+  year       = 2014
+}
+
+@InProceedings{y_li,
+      author = {Y. Li and W. Ji},
+       title = {{Thermal Analysis of Pebble-Bed Reactors Based on a Tightly Coupled Mechanical-Thermal Model}},
+  booktitle  = {{Proceedings of NURETH}},
+  year       = 2016
+}
+
+@Article{gunn1987_htc,
+   author = {D.J. Gunn},
+    title = {{Transfer of Heat or Mass to Particles in Fixed and Fluidised Beds}},
+  journal = {International Journal of Heat and Mass Transfer},
+     year = 1978,
+   volume = 21,
+    pages = {467--476},
+      DOI = {10.1016/0017-9310(78)90080-7}
+}
+
+@Article{littman,
+   author = {H. Littman and R.G. Barile and A.H. Pulsifer},
+    title = {{Gas-Particle Heat Transfer Coefficients in Packed Beds at Low {Reynolds} Numbers}},
+  journal = {I \& EC Technology},
+     year = 1968,
+   volume = 7,
+    pages = {554--561},
+      DOI = {10.1021/i160028a005}
+}
+
+@InProceedings{becker,
+   author = {S. Becker and E. Laurien},
+    title = {{Three-Dimensional Numerical Simulation of Flow and Heat Transport in High-Temperature Nuclear Reactors}},
+  booktitle  = {{High Performance Computing in Science and Engineering}},
+  year       = 2002,
+        DOI = {10.1016/S0029-5493(03)00011-6}
+}
diff --git a/modules/navier_stokes/doc/content/modules/navier_stokes/pinsfv.md b/modules/navier_stokes/doc/content/modules/navier_stokes/pinsfv.md
index 1c1247dfb9a2..a7b5b237db21 100644
--- a/modules/navier_stokes/doc/content/modules/navier_stokes/pinsfv.md
+++ b/modules/navier_stokes/doc/content/modules/navier_stokes/pinsfv.md
@@ -1,20 +1,24 @@
 # Finite Volume Incompressible Porous media Navier Stokes
 
+!alert note
+The weakly compressible formulation is also implemented for porous flow finite volume.
+
 ## Equations
 
 This module implements the porous media Navier Stokes equations. They are expressed in terms of the superficial
-velocity $\vec{v}_d = \epsilon \vec{V}$ where $\epsilon$ is the porosity and $\vec{V}$ the interstitial velocity. The
+velocity $\vec{v}_D = \epsilon \vec{v}$ where $\epsilon$ is the porosity and $\vec{v}$ the interstitial velocity. The
 superficial velocity is also known as the extrinsic or Darcy velocity. The other non-linear variables used are
 pressure and temperature. This is known as the primitive superficial set of variables.
 
 Mass equation:
 \begin{equation}
-\nabla \cdot \rho \vec{v}_d = 0
+\nabla \cdot \rho \vec{v}_D = 0
 \end{equation}
 
 Momentum equation, with friction and gravity force as example forces:
 \begin{equation}
 \dfrac{\partial \rho \mathbf{v}_D}{\partial t} + \nabla \cdot (\dfrac{\rho}{\epsilon} \mathbf{v}_D \otimes \mathbf{v}_D) = \nabla \cdot (\mu \nabla \dfrac{\mathbf{v}_D}{\epsilon}) - \epsilon \nabla p + \epsilon (\mathbf{F}_g + \mathbf{F}_f)
+\label{eq:pinsfv_mom}
 \end{equation}
 
 Fluid phase energy equation, with a convective heat transfer term:
@@ -27,7 +31,7 @@ Solid phase energy equation, with convective heat transfer and an energy source
 \dfrac{\partial (1-\epsilon) \rho c_{ps} T_s}{\partial t} = \nabla \cdot (\kappa_s \nabla T_s) + \alpha (T_f - T_s) + (1-\epsilon) \dot{Q}
 \end{equation}
 
-where $\rho$ is the fluid density, $\mu$ the viscosity, $c_p$ the specific heat capacity $\alpha$ the convective heat transfer coefficient.
+where $\rho$ is the fluid density, $\mu$ the dynamic viscosity, $c_p$ the specific heat capacity $\alpha$ the convective heat transfer coefficient. The effective diffusivities $\kappa$ include the effect of porosity, which is often approximated as $\epsilon k$ with $k$ the thermal diffusivity.
 
 ## Implementation for a straight channel
 
diff --git a/modules/navier_stokes/doc/content/source/functormaterials/FunctorErgunDragCoefficients.md b/modules/navier_stokes/doc/content/source/functormaterials/FunctorErgunDragCoefficients.md
new file mode 100644
index 000000000000..45baafc8ca14
--- /dev/null
+++ b/modules/navier_stokes/doc/content/source/functormaterials/FunctorErgunDragCoefficients.md
@@ -0,0 +1,33 @@
+# FunctorErgunDragCoefficients
+
+!syntax description /Materials/FunctorErgunDragCoefficients
+
+This class implements the Darcy and Forchheimer coefficients for the Ergun drag
+force model which is discussed in detail in
+[PINSFVMomentumFriction.md#friction_example]. We also give details on the
+definition of Darcy and Forchheimer coefficients there.
+
+To summarize, this class implements the Darcy coefficient as
+
+\begin{equation}
+\frac{150}{D_h^2}
+\end{equation}
+
+where $D_h$ is the hydraulic diameter defined as $\frac{\epsilon d_p}{1 - \epsilon}$
+where $d_p$ is the diameter of the pebbles in the bed. The Forchheimer
+coefficient is given as
+
+\begin{equation}
+\frac{2 \cdot 1.75}{D_h}
+\end{equation}
+
+where we have made the $2(1.75)$ multiplication explicit, instead of writing
+$3.5$, to make the 1.75 factor from the
+[Ergun wikipedia page](https://en.wikipedia.org/wiki/Ergun_equation) more
+recognizable.
+
+!syntax parameters /Materials/FunctorErgunDragCoefficients
+
+!syntax inputs /Materials/FunctorErgunDragCoefficients
+
+!syntax children /Materials/FunctorErgunDragCoefficients
diff --git a/modules/navier_stokes/doc/content/source/functormaterials/FunctorKappaFluid.md b/modules/navier_stokes/doc/content/source/functormaterials/FunctorKappaFluid.md
new file mode 100644
index 000000000000..4e968c373969
--- /dev/null
+++ b/modules/navier_stokes/doc/content/source/functormaterials/FunctorKappaFluid.md
@@ -0,0 +1,23 @@
+# FunctorKappaFluid
+
+!syntax description /Materials/FunctorKappaFluid
+
+## Description
+
+Most macroscale models neglect thermal dispersion [!cite](suikkanen,y_li), in which case $\kappa_f$ is given as
+
+\begin{equation}
+\label{eq-KappaFluid}
+\kappa_f=\epsilon k_f\ .
+\end{equation}
+
+Neglecting thermal dispersion is expected to be a reasonable approximation for high Reynolds numbers [!cite](gunn1987_htc,littman),
+but for low Reynolds numbers more sophisticated models should be used [!cite](becker).
+Because thermal dispersion acts to increase the diffusive effects, neglecting thermal dispersion is
+(thermally) conservative in the sense that peak temperatures are usually higher [!cite](becker).
+
+!syntax parameters /Materials/FunctorKappaFluid
+
+!syntax inputs /Materials/FunctorKappaFluid
+
+!syntax children /Materials/FunctorKappaFluid
diff --git a/modules/navier_stokes/doc/content/source/fvkernels/PINSFVMomentumFriction.md b/modules/navier_stokes/doc/content/source/fvkernels/PINSFVMomentumFriction.md
index 49fd8b560b31..8c4b03237f43 100644
--- a/modules/navier_stokes/doc/content/source/fvkernels/PINSFVMomentumFriction.md
+++ b/modules/navier_stokes/doc/content/source/fvkernels/PINSFVMomentumFriction.md
@@ -2,19 +2,115 @@
 
 This kernel adds the friction term to the porous media Navier Stokes momentum
 equations. This kernel must be used with the canonical PINSFV variable set,
-e.g. pressure and superficial velocity.  A variety of models are supported for
-the friction force:
+e.g. pressure and superficial velocity, and supports Darcy and
+Forchheimer friction models:
 
 Darcy drag model
 \begin{equation}
-F_i = - \dfrac{f_i}{\epsilon} \rho v_d
+\label{darcy}
+\epsilon F_i = - f_i \mu \epsilon \frac{v_{D,i}}{\epsilon} = -f_i \mu v_{D,i}
 \end{equation}
 Forchheimer drag model
 \begin{equation}
-F_i = - \dfrac{f_i}{\epsilon} \rho v_d
+\label{forchheimer}
+\epsilon F_i = - f_i \frac{\rho}{2} \epsilon \frac{v_{D,i}}{\epsilon}\frac{|\vec{v_D}|}{\epsilon} = - f_i \frac{\rho}{2} v_{D,i} \frac{|\vec{v_D}|}{\epsilon}
 \end{equation}
-where $F_i$ is the i-th component of the friction force, f the friction factor, which may be anisotropic,
-$\epsilon$ the porosity and $\rho$ the fluid density and $v_d$ the fluid superficial velocity.
+
+where $F_i$ is the i-th component of the friction force (denoted by
+$\mathbf{F_f}$ in [!eqref](pinsfv.md#eq:pinsfv_mom)), $f_i$ the friction factor,
+which may be anisotropic, $\epsilon$ the porosity, $\mu$ the fluid dynamic
+viscosity, $\rho$ the fluid density, and $v_{D,i}$ the i-th component of the
+fluid superficial velocity. We have used a negative sign to match the notation
+used in [!eqref](pinsfv.md#eq:pinsfv_mom) where the friction force is on the
+right-hand-side of the equation. When moved to the left-hand side, which is done
+when setting up a Newton scheme, the term becomes positive which is what is
+shown in the source code itself.  Darcy and Forchheimer terms represent
+fundamentally different friction effects. Darcy is meant to represent viscous
+effects and as shown in [darcy] has a linear dependence on the fluid
+velocity. Meanwhile, Forchheimer is meant to represent inertial effects and as
+shown in [forchheimer] has a quadratic dependence on velocity.
+
+## Computation of friction factors and pre-factors id=friction_example
+
+To outline how friction factors for Darcy and Forchheimer may be calculated,
+let's consider a specific example. We'll draw from the Ergun equation, which is
+outlined [here](https://en.wikipedia.org/wiki/Ergun_equation). Let's consider
+the form:
+
+\begin{equation}
+\Delta p = \frac{150\mu L}{d_p^2} \frac{(1-\epsilon)^2}{\epsilon^3} v_D + \frac{1.75 L \rho}{d_p} \frac{(1-\epsilon)}{\epsilon^3} |v_D| v_D
+\end{equation}
+
+where $L$ is the bed length, $\mu$ is the fluid dynamic viscosity and $d_p$ is
+representative of the diameter of the pebbles in the pebble bed. We can divide
+the equation through by $L$, recognize that $\Delta p$ denotes $p_0 - p_L$ such
+that $\Delta p/L \approx -\nabla p$, multiply the equation through by
+$-\epsilon$, move all terms to the left-hand-side, and do
+some term manipulation in order to yield:
+
+\begin{equation}
+\epsilon \nabla p + 150\mu\epsilon\frac{(1-\epsilon)^2}{\epsilon^2 d_p^2} \frac{\vec{v_D}}{\epsilon} +
+1.75\epsilon\frac{1-\epsilon}{\epsilon d_p} \rho \frac{|\vec{v}_D|}{\epsilon} \frac{\vec{v_D}}{\epsilon} = 0
+\end{equation}
+
+If we define the hydraulic diameter as $D_h = \frac{\epsilon d_p}{1 - \epsilon}$,
+then the above equation can be rewritten as:
+
+\begin{equation}
+\epsilon \nabla p + \frac{150\mu}{D_h^2} \vec{v_D} +
+\frac{1.75}{D_h} \rho \vec{v_D} \frac{|\vec{v}_D|}{\epsilon} = 0
+\end{equation}
+
+Let's introduce the interstitial fluid velocity $v = \vec{v_D} / \epsilon$ to rewrite
+the above equation as:
+
+\begin{equation}
+\epsilon \nabla p + \epsilon\frac{150\mu}{D_h^2} \vec{v} +
+\epsilon\frac{1.75}{D_h} \rho \vec{v} |\vec{v}| = 0
+\end{equation}
+
+Then dividing through by $\epsilon$:
+
+\begin{equation}
+\label{derived_ergun}
+\nabla p + \frac{150\mu}{D_h^2} \vec{v} +
+\frac{1.75}{D_h} \rho \vec{v} |\vec{v}| = 0
+\end{equation}
+
+We are now very close the forms for Darcy and Forchheimer espoused by
+[Holzmann](https://holzmann-cfd.com/community/blog-and-tools/darcy-forchheimer)
+and
+[SimScale](https://www.simscale.com/knowledge-base/predict-darcy-and-forchheimer-coefficients-for-perforated-plates-using-analytical-approach/)
+which is:
+
+\begin{equation}
+\label{holzmann}
+\nabla p + \mu D \vec{v} + \frac{\rho}{2} F |\vec{v}| \vec{v}
+\end{equation}
+
+Looking at [holzmann] we can rearrange [derived_ergun]:
+
+\begin{equation}
+\nabla p + \mu \frac{150}{D_h^2} \vec{v} +
+\frac{\rho}{2} \frac{2(1.75)}{D_h} \vec{v} |\vec{v}| = 0
+\end{equation}
+
+and arrive at the Ergun expression for the Darcy coefficient:
+
+\begin{equation}
+\frac{150}{D_h^2}
+\end{equation}
+
+and the Ergun expression for the Forchheimer coefficient:
+
+\begin{equation}
+\frac{2 \cdot 1.75}{D_h}
+\end{equation}
+
+where we have made the $2 \cdot 1.75$ multiplication explicit to make the 1.75 factor
+from the [Ergun wikipedia page](https://en.wikipedia.org/wiki/Ergun_equation)
+more recognizable. We perform a similar separation in the implementation of the
+Ergun Forchheimer coefficient outlined in [FunctorErgunDragCoefficients.md].
 
 !syntax parameters /FVKernels/PINSFVMomentumFriction
 
diff --git a/modules/navier_stokes/include/base/NS.h b/modules/navier_stokes/include/base/NS.h
index f2e0970395df..7332b4e3edce 100644
--- a/modules/navier_stokes/include/base/NS.h
+++ b/modules/navier_stokes/include/base/NS.h
@@ -22,6 +22,7 @@ using namespace HeatConduction;
 static const std::string directions[3] = {"x", "y", "z"};
 
 // geometric quantities
+static const std::string pebble_diameter = "pebble_diameter";
 static const std::string infinite_porosity = "infinite_porosity";
 static const std::string axis = "axis";
 static const std::string center = "center";
diff --git a/modules/navier_stokes/include/base/NSFVBase.h b/modules/navier_stokes/include/base/NSFVBase.h
index 1a718fdd39cb..7fadcfd3465a 100644
--- a/modules/navier_stokes/include/base/NSFVBase.h
+++ b/modules/navier_stokes/include/base/NSFVBase.h
@@ -1943,7 +1943,6 @@ NSFVBase<BaseType>::addINSMomentumFrictionKernels()
     params.template set<MooseFunctorName>(NS::density) = _density_name;
     params.template set<UserObjectName>("rhie_chow_user_object") =
         prefix() + "pins_rhie_chow_interpolator";
-    params.template set<MooseFunctorName>(NS::porosity) = _flow_porosity_functor_name;
 
     for (unsigned int block_i = 0; block_i < num_used_blocks; ++block_i)
     {
@@ -1967,10 +1966,16 @@ NSFVBase<BaseType>::addINSMomentumFrictionKernels()
         {
           const auto upper_name = MooseUtils::toUpper(_friction_types[block_i][type_i]);
           if (upper_name == "DARCY")
+          {
+            params.template set<MooseFunctorName>(NS::mu) = _dynamic_viscosity_name;
             params.template set<MooseFunctorName>("Darcy_name") = _friction_coeffs[block_i][type_i];
+          }
           else if (upper_name == "FORCHHEIMER")
+          {
             params.template set<MooseFunctorName>("Forchheimer_name") =
                 _friction_coeffs[block_i][type_i];
+            params.template set<MooseFunctorName>(NS::speed) = NS::speed;
+          }
         }
 
         getProblem().addFVKernel(kernel_type,
@@ -1990,7 +1995,6 @@ NSFVBase<BaseType>::addINSMomentumFrictionKernels()
         corr_params.template set<MooseFunctorName>(NS::density) = _density_name;
         corr_params.template set<UserObjectName>("rhie_chow_user_object") =
             prefix() + "pins_rhie_chow_interpolator";
-        corr_params.template set<MooseFunctorName>(NS::porosity) = _flow_porosity_functor_name;
         corr_params.template set<Real>("consistent_scaling") =
             parameters().template get<Real>("consistent_scaling");
         for (unsigned int d = 0; d < _dim; ++d)
@@ -2001,11 +2005,17 @@ NSFVBase<BaseType>::addINSMomentumFrictionKernels()
           {
             const auto upper_name = MooseUtils::toUpper(_friction_types[block_i][type_i]);
             if (upper_name == "DARCY")
+            {
+              corr_params.template set<MooseFunctorName>(NS::mu) = _dynamic_viscosity_name;
               corr_params.template set<MooseFunctorName>("Darcy_name") =
                   _friction_coeffs[block_i][type_i];
+            }
             else if (upper_name == "FORCHHEIMER")
+            {
               corr_params.template set<MooseFunctorName>("Forchheimer_name") =
                   _friction_coeffs[block_i][type_i];
+              corr_params.template set<MooseFunctorName>(NS::speed) = NS::speed;
+            }
           }
 
           getProblem().addFVKernel(correction_kernel_type,
@@ -2875,14 +2885,15 @@ template <class BaseType>
 void
 NSFVBase<BaseType>::addEnthalpyMaterial()
 {
-  InputParameters params = getFactory().getValidParams("INSFVEnthalpyMaterial");
+  InputParameters params = getFactory().getValidParams("INSFVEnthalpyFunctorMaterial");
   assignBlocks(params, _blocks);
 
   params.template set<MooseFunctorName>(NS::density) = _density_name;
   params.template set<MooseFunctorName>(NS::cp) = _specific_heat_name;
   params.template set<MooseFunctorName>("temperature") = _fluid_temperature_name;
 
-  getProblem().addMaterial("INSFVEnthalpyMaterial", prefix() + "ins_enthalpy_material", params);
+  getProblem().addMaterial(
+      "INSFVEnthalpyFunctorMaterial", prefix() + "ins_enthalpy_material", params);
 }
 
 template <class BaseType>
diff --git a/modules/navier_stokes/include/functormaterials/FunctorDragCoefficients.h b/modules/navier_stokes/include/functormaterials/FunctorDragCoefficients.h
new file mode 100644
index 000000000000..4534fa5d0751
--- /dev/null
+++ b/modules/navier_stokes/include/functormaterials/FunctorDragCoefficients.h
@@ -0,0 +1,25 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "FunctorMaterial.h"
+
+/**
+ * Abstract base class material providing the drag coefficients for linear and quadratic friction
+ * models. The expression of the coefficients is highly dependent on the formulation of the kernel.
+ * The reader should consult the PINSFVMomentumFrictionKernel documentation for details
+ */
+class FunctorDragCoefficients : public FunctorMaterial
+{
+public:
+  FunctorDragCoefficients(const InputParameters & parameters);
+
+  static InputParameters validParams();
+};
diff --git a/modules/navier_stokes/include/functormaterials/FunctorEffectiveFluidThermalConductivity.h b/modules/navier_stokes/include/functormaterials/FunctorEffectiveFluidThermalConductivity.h
new file mode 100644
index 000000000000..28d4a3af5bc4
--- /dev/null
+++ b/modules/navier_stokes/include/functormaterials/FunctorEffectiveFluidThermalConductivity.h
@@ -0,0 +1,29 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "FunctorMaterial.h"
+/**
+ * This is a base class material to calculate the effective thermal
+ * conductivity of the fluid phase. Unlike for the solid effective thermal
+ * conductivity, a factor of porosity is removed such that all parameters in the
+ * fluid effective thermal conductivity calculation can be expressed as porosity
+ * multiplying some parameter as $\kappa_f=\epsilon K$, where this material
+ * actually computes $K$. This is performed such that the correct spatial
+ * derivative of $\kappa_f$ can be performed, since the porosity is in general a
+ * function of space, and spatial derivatives of materials cannot be computed.
+ */
+class FunctorEffectiveFluidThermalConductivity : public FunctorMaterial
+{
+public:
+  FunctorEffectiveFluidThermalConductivity(const InputParameters & parameters);
+
+  static InputParameters validParams();
+};
diff --git a/modules/navier_stokes/include/functormaterials/FunctorErgunDragCoefficients.h b/modules/navier_stokes/include/functormaterials/FunctorErgunDragCoefficients.h
new file mode 100644
index 000000000000..122944c1bd02
--- /dev/null
+++ b/modules/navier_stokes/include/functormaterials/FunctorErgunDragCoefficients.h
@@ -0,0 +1,50 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "FunctorPebbleBedDragCoefficients.h"
+
+/**
+ *  Material providing the interphase drag coefficient according to the correlation
+ *  provided by \cite ergun.
+ *  NOTE: Do not inherit this class, inherit FunctorPebbleBedDragCoefficients or turn this
+ *       into a template before inheriting, for the sake of the CRTP static polymorphism
+ */
+class FunctorErgunDragCoefficients final
+  : public FunctorPebbleBedDragCoefficients<FunctorErgunDragCoefficients>
+{
+  friend class FunctorIsotropicDragCoefficients<FunctorErgunDragCoefficients>;
+
+public:
+  FunctorErgunDragCoefficients(const InputParameters & parameters);
+
+  static InputParameters validParams();
+
+protected:
+  template <typename Space, typename Time>
+  ADReal computeDarcyCoefficient(const Space & r, const Time & t) const;
+
+  template <typename Space, typename Time>
+  ADReal computeForchheimerCoefficient(const Space & r, const Time & t) const;
+};
+
+template <typename Space, typename Time>
+ADReal
+FunctorErgunDragCoefficients::computeDarcyCoefficient(const Space &, const Time &) const
+{
+  return 150.0;
+}
+
+template <typename Space, typename Time>
+ADReal
+FunctorErgunDragCoefficients::computeForchheimerCoefficient(const Space &, const Time &) const
+{
+  return 1.75;
+}
diff --git a/modules/navier_stokes/include/functormaterials/FunctorIsotropicDragCoefficients.h b/modules/navier_stokes/include/functormaterials/FunctorIsotropicDragCoefficients.h
new file mode 100644
index 000000000000..34a319724aee
--- /dev/null
+++ b/modules/navier_stokes/include/functormaterials/FunctorIsotropicDragCoefficients.h
@@ -0,0 +1,74 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "FunctorDragCoefficients.h"
+#include "NS.h"
+
+/**
+ * Abstract base class to compute isotropic drag coefficients, where \f$C_L\f$ and
+ * \f$C_Q\f$ are independent of direction. Each of \f$C_L\f$ and\f$C_Q\f$
+ * are separated into the calculation of a prefactor and a coefficient,where the
+ * total drag coefficient is the multiplication of these two. This separation
+ * reduces some code duplication.
+ */
+template <typename Derived>
+class FunctorIsotropicDragCoefficients : public FunctorDragCoefficients
+{
+public:
+  FunctorIsotropicDragCoefficients(const InputParameters & parameters);
+
+  static InputParameters validParams();
+
+protected:
+  /// a multiplier for adjusting Darcy coefficients
+  Real _darcy_mult;
+
+  /// a multiplier for adjusting Forchheimer coefficients
+  Real _forchheimer_mult;
+};
+
+template <typename Derived>
+InputParameters
+FunctorIsotropicDragCoefficients<Derived>::validParams()
+{
+  auto params = FunctorDragCoefficients::validParams();
+  params.addParam<Real>("Darcy_multiplier", 1, "A multiplier to adjust Darcy coefficients");
+  params.addParam<Real>(
+      "Forchheimer_multiplier", 1, "A multiplier to adjust Forchheimer coefficients");
+  return params;
+}
+
+template <typename Derived>
+FunctorIsotropicDragCoefficients<Derived>::FunctorIsotropicDragCoefficients(
+    const InputParameters & parameters)
+  : FunctorDragCoefficients(parameters),
+    _darcy_mult(getParam<Real>("Darcy_multiplier")),
+    _forchheimer_mult(getParam<Real>("Forchheimer_multiplier"))
+{
+  addFunctorProperty<ADRealVectorValue>(
+      NS::cL,
+      [this](const auto & r, const auto & t) -> ADRealVectorValue
+      {
+        RealVectorValue multipliers(1.0, 1.0, 1.0);
+        return multipliers * static_cast<Derived *>(this)->computeDarcyPrefactor(r, t) *
+               static_cast<Derived *>(this)->computeDarcyCoefficient(r, t) * _darcy_mult;
+      });
+
+  addFunctorProperty<ADRealVectorValue>(
+      NS::cQ,
+      [this](const auto & r, const auto & t) -> ADRealVectorValue
+      {
+        RealVectorValue multipliers(1.0, 1.0, 1.0);
+        return multipliers * static_cast<Derived *>(this)->computeForchheimerPrefactor(r, t) *
+               static_cast<Derived *>(this)->computeForchheimerCoefficient(r, t) *
+               _forchheimer_mult;
+      });
+}
diff --git a/modules/navier_stokes/include/functormaterials/FunctorIsotropicEffectiveFluidThermalConductivity.h b/modules/navier_stokes/include/functormaterials/FunctorIsotropicEffectiveFluidThermalConductivity.h
new file mode 100644
index 000000000000..69728ec94430
--- /dev/null
+++ b/modules/navier_stokes/include/functormaterials/FunctorIsotropicEffectiveFluidThermalConductivity.h
@@ -0,0 +1,47 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "FunctorEffectiveFluidThermalConductivity.h"
+#include "NS.h"
+
+/**
+ * Material providing an isotropic effective fluid thermal conductivity.
+ */
+template <typename Derived>
+class FunctorIsotropicEffectiveFluidThermalConductivity
+  : public FunctorEffectiveFluidThermalConductivity
+{
+public:
+  FunctorIsotropicEffectiveFluidThermalConductivity(const InputParameters & parameters);
+
+  static InputParameters validParams();
+};
+
+template <typename Derived>
+InputParameters
+FunctorIsotropicEffectiveFluidThermalConductivity<Derived>::validParams()
+{
+  return FunctorEffectiveFluidThermalConductivity::validParams();
+}
+
+template <typename Derived>
+FunctorIsotropicEffectiveFluidThermalConductivity<
+    Derived>::FunctorIsotropicEffectiveFluidThermalConductivity(const InputParameters & parameters)
+  : FunctorEffectiveFluidThermalConductivity(parameters)
+{
+  addFunctorProperty<ADRealVectorValue>(
+      NS::kappa,
+      [this](const auto & r, const auto & t) -> ADRealVectorValue
+      {
+        RealVectorValue multipliers(1.0, 1.0, 1.0);
+        return multipliers * static_cast<Derived *>(this)->computeEffectiveConductivity(r, t);
+      });
+}
diff --git a/modules/navier_stokes/include/functormaterials/FunctorKappaFluid.h b/modules/navier_stokes/include/functormaterials/FunctorKappaFluid.h
new file mode 100644
index 000000000000..c0c12e3fa108
--- /dev/null
+++ b/modules/navier_stokes/include/functormaterials/FunctorKappaFluid.h
@@ -0,0 +1,72 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "FunctorIsotropicEffectiveFluidThermalConductivity.h"
+
+/**
+ * Template for the fluid effective thermal conductivity based on a
+ * volume-average of the thermal conductivity as $\kappa_f=\epsilon k_f$. In other
+ * words, no additional effects of thermal dispersion are considered here.
+ */
+template <typename Derived>
+class FunctorKappaFluidTempl : public FunctorIsotropicEffectiveFluidThermalConductivity<Derived>
+{
+  friend class FunctorIsotropicEffectiveFluidThermalConductivity<Derived>;
+
+public:
+  FunctorKappaFluidTempl(const InputParameters & parameters);
+
+  static InputParameters validParams();
+
+protected:
+  template <typename Space, typename Time>
+  ADReal computeEffectiveConductivity(const Space & r, const Time & t) const;
+
+  /// Fluid thermal conductivity
+  const Moose::Functor<ADReal> & _k;
+};
+
+template <typename Derived>
+template <typename Space, typename Time>
+ADReal
+FunctorKappaFluidTempl<Derived>::computeEffectiveConductivity(const Space & r, const Time & t) const
+{
+  return _k(r, t);
+}
+
+/**
+ * Instantiation of the template for the fluid effective thermal conductivity based on a
+ * volume-average of the thermal conductivity as $\kappa_f=\epsilon k_f$. In other
+ * words, no additional effects of thermal dispersion are considered here.
+ * NOTE: Do not inherit this class, inherit FunctorKappaFluidTempl for the sake of the CRTP
+ *       static polymorphism
+ */
+class FunctorKappaFluid final : public FunctorKappaFluidTempl<FunctorKappaFluid>
+{
+public:
+  FunctorKappaFluid(const InputParameters & parameters);
+};
+
+template <typename Derived>
+FunctorKappaFluidTempl<Derived>::FunctorKappaFluidTempl(const InputParameters & parameters)
+  : FunctorIsotropicEffectiveFluidThermalConductivity<Derived>(parameters),
+    _k(FunctorMaterial::getFunctor<ADReal>(NS::k))
+{
+}
+
+template <typename Derived>
+InputParameters
+FunctorKappaFluidTempl<Derived>::validParams()
+{
+  auto params = FunctorIsotropicEffectiveFluidThermalConductivity<Derived>::validParams();
+  params.addClassDescription("Zero-thermal dispersion conductivity");
+  return params;
+}
diff --git a/modules/navier_stokes/include/functormaterials/FunctorPebbleBedDragCoefficients.h b/modules/navier_stokes/include/functormaterials/FunctorPebbleBedDragCoefficients.h
new file mode 100644
index 000000000000..7965b114d592
--- /dev/null
+++ b/modules/navier_stokes/include/functormaterials/FunctorPebbleBedDragCoefficients.h
@@ -0,0 +1,108 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "FunctorIsotropicDragCoefficients.h"
+#include "NS.h"
+
+/**
+ * Abstract base class to compute isotropic drag coefficients in a
+ * pebble bed. These correlations generally share common proportionalities to
+ * porosity, hydraulic diameter, etc. that are placed into this
+ * base class. The Darcy coefficient in this case scales a term with
+ * proportionality \f$\frac{1}{D_h^2}\f$,
+ * while the Forchheimer coefficient scales a term with proportionality
+ * \f$\frac{1}{D_h}\f$.
+ */
+template <typename Derived>
+class FunctorPebbleBedDragCoefficients : public FunctorIsotropicDragCoefficients<Derived>
+{
+  friend class FunctorIsotropicDragCoefficients<Derived>;
+
+public:
+  FunctorPebbleBedDragCoefficients(const InputParameters & parameters);
+
+  static InputParameters validParams();
+
+protected:
+  template <typename Space, typename Time>
+  ADReal computeDarcyPrefactor(const Space & r, const Time & t);
+
+  template <typename Space, typename Time>
+  ADReal computeForchheimerPrefactor(const Space & r, const Time & t);
+
+  /// Compute hydraulic diameter in the bed
+  template <typename Space, typename Time>
+  ADReal computeHydraulicDiameter(const Space & r, const Time & t);
+
+  /// porosity
+  const Moose::Functor<ADReal> & _eps;
+
+  /// fluid density
+  const Moose::Functor<ADReal> & _rho;
+
+  /// pebble diameter
+  const Real & _d_pebble;
+};
+
+template <typename Derived>
+InputParameters
+FunctorPebbleBedDragCoefficients<Derived>::validParams()
+{
+  InputParameters params = FunctorIsotropicDragCoefficients<Derived>::validParams();
+  params.addRequiredRangeCheckedParam<Real>(
+      NS::pebble_diameter, NS::pebble_diameter + " > 0.0", "Pebble diameter");
+  params.addRequiredParam<MooseFunctorName>(NS::porosity, "Porosity");
+
+  params.addParam<MooseFunctorName>(NS::density, NS::density, "Density");
+  return params;
+}
+
+template <typename Derived>
+FunctorPebbleBedDragCoefficients<Derived>::FunctorPebbleBedDragCoefficients(
+    const InputParameters & parameters)
+  : FunctorIsotropicDragCoefficients<Derived>(parameters),
+    _eps(this->template getFunctor<ADReal>(NS::porosity)),
+    _rho(this->template getFunctor<ADReal>(NS::density)),
+    _d_pebble(this->template getParam<Real>(NS::pebble_diameter))
+{
+}
+
+template <typename Derived>
+template <typename Space, typename Time>
+ADReal
+FunctorPebbleBedDragCoefficients<Derived>::computeDarcyPrefactor(const Space & r, const Time & t)
+{
+  const auto Dh = computeHydraulicDiameter(r, t);
+  return 1 / (Dh * Dh);
+}
+
+template <typename Derived>
+template <typename Space, typename Time>
+ADReal
+FunctorPebbleBedDragCoefficients<Derived>::computeForchheimerPrefactor(const Space & r,
+                                                                       const Time & t)
+{
+  // Why the factor of 2? Because we multiply by a factor of 1/2 in the friction kernel to be
+  // consistent with the formulation of Forchheimer in
+  // https://holzmann-cfd.com/community/blog-and-tools/darcy-forchheimer and
+  // https://www.simscale.com/knowledge-base/predict-darcy-and-forchheimer-coefficients-for-perforated-plates-using-analytical-approach/
+  return 2 / computeHydraulicDiameter(r, t);
+}
+
+template <typename Derived>
+template <typename Space, typename Time>
+ADReal
+FunctorPebbleBedDragCoefficients<Derived>::computeHydraulicDiameter(const Space & r, const Time & t)
+{
+  mooseAssert(1.0 - _eps(r, t).value() > 1e-8,
+              "Bed hydraulic diameter is ill-defined at porosity of 1.");
+  return _eps(r, t) * _d_pebble / (1.0 - _eps(r, t));
+}
diff --git a/modules/navier_stokes/include/fvkernels/PINSFVMomentumFriction.h b/modules/navier_stokes/include/fvkernels/PINSFVMomentumFriction.h
index fd109c21e115..67832d654bf6 100644
--- a/modules/navier_stokes/include/fvkernels/PINSFVMomentumFriction.h
+++ b/modules/navier_stokes/include/fvkernels/PINSFVMomentumFriction.h
@@ -31,14 +31,16 @@ class PINSFVMomentumFriction : public INSFVElementalKernel
                                      const Moose::StateArg & state);
 
   /// Darcy coefficient
-  const Moose::Functor<ADRealVectorValue> * const _cL;
+  const Moose::Functor<ADRealVectorValue> * const _D;
   /// Forchheimer coefficient
-  const Moose::Functor<ADRealVectorValue> * const _cQ;
+  const Moose::Functor<ADRealVectorValue> * const _F;
   /// Booleans to select the right models
   const bool _use_Darcy_friction_model;
   const bool _use_Forchheimer_friction_model;
-  /// Porosity to compute the intersitial velocity from the superficial velocity
-  const Moose::Functor<ADReal> & _eps;
+  /// Dynamic viscosity
+  const Moose::Functor<ADReal> * const _mu;
   /// Density as a functor
   const Moose::Functor<ADReal> & _rho;
+  /// Speed (norm of the interstitial velocity) as a functor
+  const Moose::Functor<ADReal> * const _speed;
 };
diff --git a/modules/navier_stokes/include/fvkernels/PINSFVMomentumFrictionCorrection.h b/modules/navier_stokes/include/fvkernels/PINSFVMomentumFrictionCorrection.h
index b43edd32de06..d103f9d1994a 100644
--- a/modules/navier_stokes/include/fvkernels/PINSFVMomentumFrictionCorrection.h
+++ b/modules/navier_stokes/include/fvkernels/PINSFVMomentumFrictionCorrection.h
@@ -21,18 +21,20 @@ class PINSFVMomentumFrictionCorrection : public INSFVFluxKernel
 
 protected:
   /// Darcy coefficient
-  const Moose::Functor<ADRealVectorValue> * const _cL;
+  const Moose::Functor<ADRealVectorValue> * const _D;
   /// Forchheimer coefficient
-  const Moose::Functor<ADRealVectorValue> * const _cQ;
+  const Moose::Functor<ADRealVectorValue> * const _F;
 
   /// Booleans to select the right models
   const bool _use_Darcy_friction_model;
   const bool _use_Forchheimer_friction_model;
 
-  /// Porosity to compute the interstitial velocity from the superficial velocity
-  const Moose::Functor<ADReal> & _eps;
   /// Density as a functor
   const Moose::Functor<ADReal> & _rho;
+  /// Dynamic viscosity as a functor
+  const Moose::Functor<ADReal> & _mu;
+  /// Speed (norm of the interstitial velocity) as a functor
+  const Moose::Functor<ADReal> * const _speed;
 
   /// Parameter for scaling the consistent pressure interpolation
   Real _consistent_scaling;
diff --git a/modules/navier_stokes/src/functormaterials/FunctorDragCoefficients.C b/modules/navier_stokes/src/functormaterials/FunctorDragCoefficients.C
new file mode 100644
index 000000000000..d7317859bb90
--- /dev/null
+++ b/modules/navier_stokes/src/functormaterials/FunctorDragCoefficients.C
@@ -0,0 +1,21 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "FunctorDragCoefficients.h"
+
+InputParameters
+FunctorDragCoefficients::validParams()
+{
+  return FunctorMaterial::validParams();
+}
+
+FunctorDragCoefficients::FunctorDragCoefficients(const InputParameters & parameters)
+  : FunctorMaterial(parameters)
+{
+}
diff --git a/modules/navier_stokes/src/functormaterials/FunctorEffectiveFluidThermalConductivity.C b/modules/navier_stokes/src/functormaterials/FunctorEffectiveFluidThermalConductivity.C
new file mode 100644
index 000000000000..6ec031629a5c
--- /dev/null
+++ b/modules/navier_stokes/src/functormaterials/FunctorEffectiveFluidThermalConductivity.C
@@ -0,0 +1,24 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "FunctorEffectiveFluidThermalConductivity.h"
+
+InputParameters
+FunctorEffectiveFluidThermalConductivity::validParams()
+{
+  auto params = FunctorMaterial::validParams();
+  params.addClassDescription("Base class for computing effective fluid thermal conductivity");
+  return params;
+}
+
+FunctorEffectiveFluidThermalConductivity::FunctorEffectiveFluidThermalConductivity(
+    const InputParameters & parameters)
+  : FunctorMaterial(parameters)
+{
+}
diff --git a/modules/navier_stokes/src/functormaterials/FunctorErgunDragCoefficients.C b/modules/navier_stokes/src/functormaterials/FunctorErgunDragCoefficients.C
new file mode 100644
index 000000000000..4a483fece607
--- /dev/null
+++ b/modules/navier_stokes/src/functormaterials/FunctorErgunDragCoefficients.C
@@ -0,0 +1,26 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "FunctorErgunDragCoefficients.h"
+
+registerMooseObject("NavierStokesApp", FunctorErgunDragCoefficients);
+
+InputParameters
+FunctorErgunDragCoefficients::validParams()
+{
+  auto params = FunctorPebbleBedDragCoefficients<FunctorErgunDragCoefficients>::validParams();
+  params.addClassDescription("Material providing linear and quadratic drag "
+                             "coefficients based on the correlation developed by Ergun.");
+  return params;
+}
+
+FunctorErgunDragCoefficients::FunctorErgunDragCoefficients(const InputParameters & parameters)
+  : FunctorPebbleBedDragCoefficients<FunctorErgunDragCoefficients>(parameters)
+{
+}
diff --git a/modules/navier_stokes/src/functormaterials/FunctorKappaFluid.C b/modules/navier_stokes/src/functormaterials/FunctorKappaFluid.C
new file mode 100644
index 000000000000..fe6fbadb7c09
--- /dev/null
+++ b/modules/navier_stokes/src/functormaterials/FunctorKappaFluid.C
@@ -0,0 +1,17 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "FunctorKappaFluid.h"
+
+registerMooseObject("NavierStokesApp", FunctorKappaFluid);
+
+FunctorKappaFluid::FunctorKappaFluid(const InputParameters & parameters)
+  : FunctorKappaFluidTempl<FunctorKappaFluid>(parameters)
+{
+}
diff --git a/modules/navier_stokes/src/functormaterials/PINSFVSpeedFunctorMaterial.C b/modules/navier_stokes/src/functormaterials/PINSFVSpeedFunctorMaterial.C
index aeef39f1cf51..d873fc25d374 100644
--- a/modules/navier_stokes/src/functormaterials/PINSFVSpeedFunctorMaterial.C
+++ b/modules/navier_stokes/src/functormaterials/PINSFVSpeedFunctorMaterial.C
@@ -9,6 +9,7 @@
 
 #include "PINSFVSpeedFunctorMaterial.h"
 #include "NS.h"
+#include "NavierStokesMethods.h"
 
 registerMooseObject("NavierStokesApp", PINSFVSpeedFunctorMaterial);
 
@@ -27,7 +28,18 @@ PINSFVSpeedFunctorMaterial::validParams()
       NS::superficial_velocity_y, 0, "The y component of the fluid superficial velocity variable.");
   params.addParam<MooseFunctorName>(
       NS::superficial_velocity_z, 0, "The z component of the fluid superficial velocity variable.");
-  params.addParam<MooseFunctorName>(NS::T_fluid, "The fluid temperature variable.");
+  auto add_property = [&params](const auto & property_name)
+  {
+    params.addParam<MooseFunctorName>(property_name,
+                                      property_name,
+                                      "The name to give the declared '" + property_name +
+                                          "' functor property");
+  };
+  add_property(NS::velocity);
+  add_property(NS::speed);
+  add_property(NS::velocity_x);
+  add_property(NS::velocity_y);
+  add_property(NS::velocity_z);
 
   return params;
 }
@@ -51,33 +63,28 @@ PINSFVSpeedFunctorMaterial::PINSFVSpeedFunctorMaterial(const InputParameters & p
                blocksMaxDimension());
 
   // Interstitial velocity is needed by certain correlations
-  addFunctorProperty<ADRealVectorValue>(NS::velocity,
-                                        [this](const auto & r, const auto & t) -> ADRealVectorValue
-                                        {
-                                          return ADRealVectorValue(_superficial_vel_x(r, t),
-                                                                   _superficial_vel_y(r, t),
-                                                                   _superficial_vel_z(r, t)) /
-                                                 _eps(r, t);
-                                        });
+  const auto & interstitial_velocity = addFunctorProperty<ADRealVectorValue>(
+      getParam<MooseFunctorName>(NS::velocity),
+      [this](const auto & r, const auto & t) -> ADRealVectorValue
+      {
+        return ADRealVectorValue(
+                   _superficial_vel_x(r, t), _superficial_vel_y(r, t), _superficial_vel_z(r, t)) /
+               _eps(r, t);
+      });
 
   // Speed is normal of regular interstitial velocity
   // This is needed to compute the Reynolds number
-  addFunctorProperty<ADReal>(NS::speed,
-                             [this](const auto & r, const auto & t) -> ADReal
-                             {
-                               // if the velocity is zero, then the norm function call fails because
-                               // AD tries to calculate the derivatives which causes a divide by
-                               // zero - because d/dx(sqrt(f(x))) = 1/2/sqrt(f(x))*df/dx. So add a
-                               // bit of noise to avoid this failure mode.
-                               if ((MooseUtils::absoluteFuzzyEqual(_superficial_vel_x(r, t), 0)) &&
-                                   (MooseUtils::absoluteFuzzyEqual(_superficial_vel_y(r, t), 0)) &&
-                                   (MooseUtils::absoluteFuzzyEqual(_superficial_vel_z(r, t), 0)))
-                                 return 1e-42;
+  addFunctorProperty<ADReal>(getParam<MooseFunctorName>(NS::speed),
+                             [&interstitial_velocity](const auto & r, const auto & t) -> ADReal
+                             { return NS::computeSpeed(interstitial_velocity(r, t)); });
 
-                               return ADRealVectorValue(_superficial_vel_x(r, t),
-                                                        _superficial_vel_y(r, t),
-                                                        _superficial_vel_z(r, t))
-                                          .norm() /
-                                      _eps(r, t);
-                             });
+  addFunctorProperty<ADReal>(getParam<MooseFunctorName>(NS::velocity_x),
+                             [&interstitial_velocity](const auto & r, const auto & t) -> ADReal
+                             { return interstitial_velocity(r, t)(0); });
+  addFunctorProperty<ADReal>(getParam<MooseFunctorName>(NS::velocity_y),
+                             [&interstitial_velocity](const auto & r, const auto & t) -> ADReal
+                             { return interstitial_velocity(r, t)(1); });
+  addFunctorProperty<ADReal>(getParam<MooseFunctorName>(NS::velocity_z),
+                             [&interstitial_velocity](const auto & r, const auto & t) -> ADReal
+                             { return interstitial_velocity(r, t)(2); });
 }
diff --git a/modules/navier_stokes/src/fvkernels/PINSFVMomentumFriction.C b/modules/navier_stokes/src/fvkernels/PINSFVMomentumFriction.C
index b6329ece0dad..b1f6eaa0eff1 100644
--- a/modules/navier_stokes/src/fvkernels/PINSFVMomentumFriction.C
+++ b/modules/navier_stokes/src/fvkernels/PINSFVMomentumFriction.C
@@ -24,34 +24,74 @@ PINSFVMomentumFriction::validParams()
   params.addParam<MooseFunctorName>("Darcy_name", "Name of the Darcy coefficients property.");
   params.addParam<MooseFunctorName>("Forchheimer_name",
                                     "Name of the Forchheimer coefficients property.");
-  params.addParam<MooseFunctorName>(NS::porosity, NS::porosity, "the porosity");
+  params.addParam<MooseFunctorName>(NS::mu, "The dynamic viscosity");
+  params.addParam<MooseFunctorName>(
+      NS::speed,
+      "The norm of the interstitial velocity. This is required for Forchheimer calculations");
   params.addRequiredParam<MooseFunctorName>(NS::density, "The density.");
   return params;
 }
 
 PINSFVMomentumFriction::PINSFVMomentumFriction(const InputParameters & params)
   : INSFVElementalKernel(params),
-    _cL(isParamValid("Darcy_name") ? &getFunctor<ADRealVectorValue>("Darcy_name") : nullptr),
-    _cQ(isParamValid("Forchheimer_name") ? &getFunctor<ADRealVectorValue>("Forchheimer_name")
-                                         : nullptr),
+    _D(isParamValid("Darcy_name") ? &getFunctor<ADRealVectorValue>("Darcy_name") : nullptr),
+    _F(isParamValid("Forchheimer_name") ? &getFunctor<ADRealVectorValue>("Forchheimer_name")
+                                        : nullptr),
     _use_Darcy_friction_model(isParamValid("Darcy_name")),
     _use_Forchheimer_friction_model(isParamValid("Forchheimer_name")),
-    _eps(getFunctor<ADReal>(NS::porosity)),
-    _rho(getFunctor<ADReal>(NS::density))
+    _mu(isParamValid(NS::mu) ? &getFunctor<ADReal>(NS::mu) : nullptr),
+    _rho(getFunctor<ADReal>(NS::density)),
+    _speed(isParamValid(NS::speed) ? &getFunctor<ADReal>(NS::speed) : nullptr)
 {
   if (!_use_Darcy_friction_model && !_use_Forchheimer_friction_model)
     mooseError("At least one friction model needs to be specified.");
+
+  if (_use_Forchheimer_friction_model && !_speed)
+    mooseError("If using a Forchheimer friction model, then the '",
+               NS::speed,
+               "' parameter must be provided");
+
+  if (_use_Darcy_friction_model && !_mu)
+    mooseError(
+        "If using a Darcy friction model, then the '", NS::mu, "' parameter must be provided");
 }
 
+// We are going to follow the formulations in
+// https://holzmann-cfd.com/community/blog-and-tools/darcy-forchheimer and
+// https://www.simscale.com/knowledge-base/predict-darcy-and-forchheimer-coefficients-for-perforated-plates-using-analytical-approach/
+// for Darcy and Forchheimer which define:
+//
+// \nabla p = \mu D v + \frac{\rho}{2} F |v| v
+//
+// where v denotes the interstitial velocity (e.g. the true fluid velocity). Note that to be
+// consistent with our overall porous Navier-Stokes formulation, we must multiply the above equation
+// by porosity:
+//
+// \epsilon \nabla p = \epsilon \mu D v + \epsilon \frac{\rho}{2} F |v| v =
+// \mu D v_D + \frac{\rho}{2} F |v| v_D
+//
+// where v_D is the superficial velocity, which is the variable this kernel is acting on. The two
+// terms on the final RHS are what this kernel implements. Because both monolithic and segregated
+// solves multiply by v_D in the gatherRCData and computeSegregatedContribution methods
+// respectively, it is the job of the computeFrictionWCoefficient method to compute, for the Darcy
+// term:
+//
+// \mu D
+//
+// and for the Forchheimer term
+//
+// \frac{\rho}{2} F |v|
+
 ADReal
 PINSFVMomentumFriction::computeFrictionWCoefficient(const Moose::ElemArg & elem_arg,
                                                     const Moose::StateArg & state)
 {
   ADReal coefficient = 0.0;
   if (_use_Darcy_friction_model)
-    coefficient += (*_cL)(elem_arg, state)(_index)*_rho(elem_arg, state) / _eps(elem_arg, state);
+    coefficient += (*_mu)(elem_arg, state) * (*_D)(elem_arg, state)(_index);
   if (_use_Forchheimer_friction_model)
-    coefficient += (*_cQ)(elem_arg, state)(_index)*_rho(elem_arg, state) / _eps(elem_arg, state);
+    coefficient +=
+        _rho(elem_arg, state) / 2 * (*_F)(elem_arg, state)(_index) * (*_speed)(elem_arg, state);
 
   return coefficient;
 }
diff --git a/modules/navier_stokes/src/fvkernels/PINSFVMomentumFrictionCorrection.C b/modules/navier_stokes/src/fvkernels/PINSFVMomentumFrictionCorrection.C
index c3360c7cc86b..724ab1458c13 100644
--- a/modules/navier_stokes/src/fvkernels/PINSFVMomentumFrictionCorrection.C
+++ b/modules/navier_stokes/src/fvkernels/PINSFVMomentumFrictionCorrection.C
@@ -24,8 +24,11 @@ PINSFVMomentumFrictionCorrection::validParams()
   params.addParam<MooseFunctorName>("Darcy_name", "Name of the Darcy coefficients property.");
   params.addParam<MooseFunctorName>("Forchheimer_name",
                                     "Name of the Forchheimer coefficients property.");
-  params.addParam<MooseFunctorName>(NS::porosity, NS::porosity, "The porosity");
   params.addRequiredParam<MooseFunctorName>(NS::density, "The density.");
+  params.addParam<MooseFunctorName>(
+      NS::speed,
+      "The norm of the interstitial velocity. This is required for Forchheimer calculations");
+  params.addParam<MooseFunctorName>(NS::mu, NS::mu, "The dynamic viscosity");
   params.addRangeCheckedParam<Real>("consistent_scaling",
                                     1,
                                     "consistent_scaling >= 0",
@@ -36,17 +39,23 @@ PINSFVMomentumFrictionCorrection::validParams()
 
 PINSFVMomentumFrictionCorrection::PINSFVMomentumFrictionCorrection(const InputParameters & params)
   : INSFVFluxKernel(params),
-    _cL(isParamValid("Darcy_name") ? &getFunctor<ADRealVectorValue>("Darcy_name") : nullptr),
-    _cQ(isParamValid("Forchheimer_name") ? &getFunctor<ADRealVectorValue>("Forchheimer_name")
-                                         : nullptr),
+    _D(isParamValid("Darcy_name") ? &getFunctor<ADRealVectorValue>("Darcy_name") : nullptr),
+    _F(isParamValid("Forchheimer_name") ? &getFunctor<ADRealVectorValue>("Forchheimer_name")
+                                        : nullptr),
     _use_Darcy_friction_model(isParamValid("Darcy_name")),
     _use_Forchheimer_friction_model(isParamValid("Forchheimer_name")),
-    _eps(getFunctor<ADReal>(NS::porosity)),
     _rho(getFunctor<ADReal>(NS::density)),
+    _mu(getFunctor<ADReal>(NS::mu)),
+    _speed(isParamValid(NS::speed) ? &getFunctor<ADReal>(NS::speed) : nullptr),
     _consistent_scaling(getParam<Real>("consistent_scaling"))
 {
   if (!_use_Darcy_friction_model && !_use_Forchheimer_friction_model)
     mooseError("At least one friction model needs to be specified.");
+
+  if (_use_Forchheimer_friction_model && !_speed)
+    mooseError("If using a Forchheimer friction model, then the '",
+               NS::speed,
+               "' parameter must be provided");
 }
 
 void
@@ -79,17 +88,15 @@ PINSFVMomentumFrictionCorrection::gatherRCData(const FaceInfo & fi)
   {
     if (_use_Darcy_friction_model)
     {
-      friction_term_elem +=
-          (*_cL)(elem_face, state)(_index)*_rho(elem_face, state) / _eps(elem_face, state);
-      friction_term_neighbor += (*_cL)(neighbor_face, state)(_index)*_rho(neighbor_face, state) /
-                                _eps(neighbor_face, state);
+      friction_term_elem += (*_D)(elem_face, state)(_index)*_mu(elem_face, state);
+      friction_term_neighbor += (*_D)(neighbor_face, state)(_index)*_mu(neighbor_face, state);
     }
     if (_use_Forchheimer_friction_model)
     {
       friction_term_elem +=
-          (*_cQ)(elem_face, state)(_index)*_rho(elem_face, state) / _eps(elem_face, state);
-      friction_term_neighbor += (*_cQ)(neighbor_face, state)(_index)*_rho(neighbor_face, state) /
-                                _eps(neighbor_face, state);
+          (*_F)(elem_face, state)(_index)*_rho(elem_face, state) / 2 * (*_speed)(elem_face, state);
+      friction_term_neighbor += (*_F)(neighbor_face, state)(_index)*_rho(neighbor_face, state) / 2 *
+                                (*_speed)(neighbor_face, state);
     }
 
     Point _neighbor_centroid = _face_info->neighborCentroid();
@@ -112,9 +119,10 @@ PINSFVMomentumFrictionCorrection::gatherRCData(const FaceInfo & fi)
     const auto face =
         makeFace(*_face_info, Moose::FV::limiterType(Moose::FV::InterpMethod::Average), true);
     if (_use_Darcy_friction_model)
-      friction_term_elem += (*_cL)(face, state)(_index)*_rho(face, state) / _eps(face, state);
+      friction_term_elem += (*_D)(face, state)(_index)*_mu(face, state);
     if (_use_Forchheimer_friction_model)
-      friction_term_elem += (*_cQ)(face, state)(_index)*_rho(face, state) / _eps(face, state);
+      friction_term_elem +=
+          (*_F)(face, state)(_index)*_rho(face, state) / 2 * (*_speed)(face, state);
 
     Real geometric_factor =
         _consistent_scaling * std::pow((_elem_centroid - _face_centroid).norm(), 2);
diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/friction/2d-rc-friction-action.i b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/friction/2d-rc-friction-action.i
index 367d4b348de8..f3fda2415792 100644
--- a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/friction/2d-rc-friction-action.i
+++ b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/friction/2d-rc-friction-action.i
@@ -52,6 +52,8 @@ rho = 1.1
     superficial_vel_x = vel_x
     superficial_vel_y = vel_y
     porosity = 1
+    vel_x = vel_x_mat
+    vel_y = vel_y_mat
   []
   [Re_material]
     type = ReynoldsNumberFunctorMaterial
diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/friction/2d-rc-friction.i b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/friction/2d-rc-friction.i
index 1e8ec8895462..469d83b842e1 100644
--- a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/friction/2d-rc-friction.i
+++ b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/friction/2d-rc-friction.i
@@ -165,6 +165,8 @@ velocity_interp_method = 'rc'
     superficial_vel_x = vel_x
     superficial_vel_y = vel_y
     porosity = 1
+    vel_x = vel_x_mat
+    vel_y = vel_y_mat
   []
   [Re_material]
     type = ReynoldsNumberFunctorMaterial
diff --git a/modules/navier_stokes/test/tests/finite_volume/materials/ergun/ergun.i b/modules/navier_stokes/test/tests/finite_volume/materials/ergun/ergun.i
new file mode 100644
index 000000000000..3a600b64cce1
--- /dev/null
+++ b/modules/navier_stokes/test/tests/finite_volume/materials/ergun/ergun.i
@@ -0,0 +1,428 @@
+# This file simulates flow of fluid in a porous elbow for the purpose of verifying
+# correct implementation of the various different solution variable sets. This input
+# tests correct implementation of the primitive superficial variable set. Flow enters on the top
+# and exits on the right. Because the purpose is only to test the equivalence of
+# different equation sets, no solid energy equation is included.
+
+porosity_left = 0.4
+porosity_right = 0.6
+pebble_diameter = 0.06
+
+mu = 1.81e-5 # This has been increased to avoid refining the mesh
+M = 28.97e-3
+R = 8.3144598
+
+# inlet mass flowrate, kg/s
+mdot = -10.0
+
+# inlet mass flux (superficial)
+mflux_in_superficial = ${fparse mdot / (pi * 0.5 * 0.5)}
+
+# inlet mass flux (interstitial)
+mflux_in_interstitial = ${fparse mflux_in_superficial / porosity_left}
+
+p_initial = 201325.0
+T_initial = 300.0
+rho_initial = ${fparse p_initial / T_initial * M / R}
+
+vel_y_initial = ${fparse mflux_in_interstitial / rho_initial}
+vel_x_initial = 0.0
+
+superficial_vel_y_initial = ${fparse mflux_in_superficial / rho_initial}
+superficial_vel_x_initial = 1e-12
+
+# Computation parameters
+velocity_interp_method = 'rc'
+advected_interp_method = 'upwind'
+
+# ==============================================================================
+# GEOMETRY AND MESH
+# ==============================================================================
+
+[Mesh]
+  [fmg]
+    type = FileMeshGenerator
+    file = 'ergun_in.e'
+  []
+  coord_type = RZ
+[]
+
+[UserObjects]
+  [rc]
+    type = PINSFVRhieChowInterpolator
+    u = superficial_vel_x
+    v = superficial_vel_y
+    pressure = pressure
+    porosity = porosity
+  []
+[]
+
+[GlobalParams]
+  porosity = porosity
+  pebble_diameter = ${pebble_diameter}
+  fp = fp
+
+  # rho for the kernels. Must match fluid property!
+  rho = ${rho_initial}
+
+  fv = true
+  velocity_interp_method = ${velocity_interp_method}
+  advected_interp_method = ${advected_interp_method}
+
+  # behavior at time of test creation
+  two_term_boundary_expansion = false
+  rhie_chow_user_object = 'rc'
+[]
+
+# ==============================================================================
+# VARIABLES AND KERNELS
+# ==============================================================================
+
+[Variables]
+  [pressure]
+    type = INSFVPressureVariable
+    initial_condition = ${p_initial}
+  []
+  [superficial_vel_x]
+    type = PINSFVSuperficialVelocityVariable
+    initial_condition = ${superficial_vel_x_initial}
+  []
+  [superficial_vel_y]
+    type = PINSFVSuperficialVelocityVariable
+    initial_condition = ${superficial_vel_y_initial}
+  []
+[]
+
+[FVKernels]
+  # Mass Equation.
+  [mass]
+    type = PINSFVMassAdvection
+    variable = 'pressure'
+  []
+
+  # Momentum x component equation.
+  [vel_x_time]
+    type = PINSFVMomentumTimeDerivative
+    variable = 'superficial_vel_x'
+    momentum_component = 'x'
+  []
+  [vel_x_advection]
+    type = PINSFVMomentumAdvection
+    variable = 'superficial_vel_x'
+    momentum_component = 'x'
+  []
+  [vel_x_viscosity]
+    type = PINSFVMomentumDiffusion
+    variable = 'superficial_vel_x'
+    momentum_component = 'x'
+    mu = 'mu'
+  []
+  [u_pressure]
+    type = PINSFVMomentumPressure
+    variable = 'superficial_vel_x'
+    pressure = pressure
+    momentum_component = 'x'
+  []
+  [u_friction]
+    type = PINSFVMomentumFriction
+    variable = 'superficial_vel_x'
+    Darcy_name = 'Darcy_coefficient'
+    Forchheimer_name = 'Forchheimer_coefficient'
+    momentum_component = 'x'
+    speed = speed
+    mu = 'mu'
+  []
+
+  # Momentum y component equation.
+  [vel_y_time]
+    type = PINSFVMomentumTimeDerivative
+    variable = 'superficial_vel_y'
+    momentum_component = 'y'
+  []
+  [vel_y_advection]
+    type = PINSFVMomentumAdvection
+    variable = 'superficial_vel_y'
+    momentum_component = 'y'
+  []
+  [vel_y_viscosity]
+    type = PINSFVMomentumDiffusion
+    variable = 'superficial_vel_y'
+    momentum_component = 'y'
+    mu = 'mu'
+  []
+  [v_pressure]
+    type = PINSFVMomentumPressure
+    variable = 'superficial_vel_y'
+    pressure = pressure
+    momentum_component = 'y'
+  []
+  [v_friction]
+    type = PINSFVMomentumFriction
+    variable = 'superficial_vel_y'
+    Darcy_name = 'Darcy_coefficient'
+    Forchheimer_name = 'Forchheimer_coefficient'
+    momentum_component = 'y'
+    mu = 'mu'
+    speed = speed
+  []
+  [gravity]
+    type = PINSFVMomentumGravity
+    variable = 'superficial_vel_y'
+    gravity = '0 -9.81 0'
+    momentum_component = 'y'
+  []
+[]
+
+# ==============================================================================
+# AUXVARIABLES AND AUXKERNELS
+# ==============================================================================
+
+[AuxVariables]
+  [T_fluid]
+    initial_condition = ${T_initial}
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [vel_x]
+    initial_condition = ${fparse vel_x_initial}
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [vel_y]
+    initial_condition = ${fparse vel_y_initial}
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [porosity_out]
+    type = MooseVariableFVReal
+  []
+[]
+
+[AuxKernels]
+  [vel_x]
+    type = FunctorAux
+    variable = vel_x
+    functor = vel_x_mat
+  []
+  [vel_y]
+    type = FunctorAux
+    variable = vel_y
+    functor = vel_y_mat
+  []
+  [porosity_out]
+    type = FunctorAux
+    variable = porosity_out
+    functor = porosity
+  []
+[]
+
+# ==============================================================================
+# FLUID PROPERTIES, MATERIALS AND USER OBJECTS
+# ==============================================================================
+
+[FluidProperties]
+  [fp]
+    type = IdealGasFluidProperties
+    k = 0.0
+    mu = ${mu}
+    gamma = 1.4
+    molar_mass = ${M}
+  []
+[]
+
+[FunctorMaterials]
+  [enthalpy]
+    type = INSFVEnthalpyMaterial
+    temperature = 'T_fluid'
+  []
+  [speed]
+    type = PINSFVSpeedFunctorMaterial
+    superficial_vel_x = 'superficial_vel_x'
+    superficial_vel_y = 'superficial_vel_y'
+    porosity = porosity
+    vel_x = vel_x_mat
+    vel_y = vel_y_mat
+  []
+  [kappa]
+    type = FunctorKappaFluid
+  []
+  [const_Fdrags_mat]
+    type = FunctorErgunDragCoefficients
+    porosity = porosity
+  []
+  [fluidprops]
+    type = GeneralFunctorFluidProps
+    mu_rampdown = mu_func
+    porosity = porosity
+    characteristic_length = ${pebble_diameter}
+    T_fluid = 'T_fluid'
+    pressure = 'pressure'
+    speed = 'speed'
+  []
+[]
+
+d = 0.05
+
+[Functions]
+  [mu_func]
+    type = PiecewiseLinear
+    x = '1 3 5 10 15 20'
+    y = '1e5 1e4 1e3 1e2 1e1 1'
+  []
+  [real_porosity_function]
+    type = ParsedFunction
+    expression = 'if (x < 0.6 - ${d}, ${porosity_left}, if (x > 0.6 + ${d}, ${porosity_right},
+        (x-(0.6-${d}))/(2*${d})*(${porosity_right}-${porosity_left}) + ${porosity_left}))'
+  []
+  [porosity]
+    type = ParsedFunction
+    expression = 'if (x < 0.6 - ${d}, ${porosity_left}, if (x > 0.6 + ${d}, ${porosity_right},
+        (x-(0.6-${d}))/(2*${d})*(${porosity_right}-${porosity_left}) + ${porosity_left}))'
+  []
+[]
+
+# ==============================================================================
+# BOUNDARY CONDITIONS
+# ==============================================================================
+
+[FVBCs]
+  [outlet_p]
+    type = INSFVOutletPressureBC
+    variable = 'pressure'
+    function = ${p_initial}
+    boundary = 'right'
+  []
+
+  ## No or Free slip BC
+  [free-slip-wall-x]
+    type = INSFVNaturalFreeSlipBC
+    boundary = 'bottom wall_1 wall_2 left'
+    variable = superficial_vel_x
+    momentum_component = 'x'
+  []
+  [free-slip-wall-y]
+    type = INSFVNaturalFreeSlipBC
+    boundary = 'bottom wall_1 wall_2 left'
+    variable = superficial_vel_y
+    momentum_component = 'y'
+  []
+
+  ## Symmetry
+  [symmetry-x]
+    type = PINSFVSymmetryVelocityBC
+    boundary = 'left'
+    variable = superficial_vel_x
+    u = superficial_vel_x
+    v = superficial_vel_y
+    mu = 'mu'
+    momentum_component = 'x'
+  []
+  [symmetry-y]
+    type = PINSFVSymmetryVelocityBC
+    boundary = 'left'
+    variable = superficial_vel_y
+    u = superficial_vel_x
+    v = superficial_vel_y
+    mu = 'mu'
+    momentum_component = 'y'
+  []
+  [symmetry-p]
+    type = INSFVSymmetryPressureBC
+    boundary = 'left'
+    variable = 'pressure'
+  []
+
+  ## inlet
+  [inlet_vel_x]
+    type = INSFVInletVelocityBC
+    variable = 'superficial_vel_x'
+    function = ${superficial_vel_x_initial}
+    boundary = 'top'
+  []
+  [inlet_vel_y]
+    type = INSFVInletVelocityBC
+    variable = 'superficial_vel_y'
+    function = ${superficial_vel_y_initial}
+    boundary = 'top'
+  []
+[]
+# ==============================================================================
+# EXECUTION PARAMETERS
+# ==============================================================================
+
+[Executioner]
+  type = Transient
+
+  solve_type = 'NEWTON'
+  petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -ksp_gmres_restart'
+  petsc_options_value = 'asm      lu           NONZERO                   200'
+  line_search = 'none'
+
+  # Problem time parameters
+  dtmin = 0.01
+  dtmax = 2000
+  end_time = 3000
+  # must be the same as the fluid
+
+  # Iterations parameters
+  l_max_its = 50
+  l_tol     = 1e-8
+  nl_max_its = 25
+  # nl_rel_tol = 5e-7
+  nl_abs_tol = 2e-7
+
+  # Automatic scaling
+  automatic_scaling = true
+  verbose = true
+
+  [TimeStepper]
+    type = IterationAdaptiveDT
+    dt                 = 0.025
+    cutback_factor     = 0.5
+    growth_factor      = 2.0
+  []
+
+  # Steady state detection.
+  steady_state_detection = true
+  steady_state_tolerance = 1e-7
+  steady_state_start_time = 400
+[]
+
+# ==============================================================================
+# POSTPROCESSORS DEBUG AND OUTPUTS
+# ==============================================================================
+
+[Postprocessors]
+  [mass_flow_in]
+    type = VolumetricFlowRate
+    boundary = 'top'
+    vel_x = 'superficial_vel_x'
+    vel_y = 'superficial_vel_y'
+    advected_quantity = ${rho_initial}
+    execute_on = 'INITIAL TIMESTEP_END'
+  []
+  [mass_flow_out]
+    type = VolumetricFlowRate
+    boundary = 'right'
+    vel_x = 'superficial_vel_x'
+    vel_y = 'superficial_vel_y'
+    advected_quantity = ${rho_initial}
+    execute_on = 'INITIAL TIMESTEP_END'
+  []
+  [p_in]
+    type = SideAverageValue
+    variable = pressure
+    boundary = 'top'
+  []
+  [dP]
+    type = LinearCombinationPostprocessor
+    pp_names = 'p_in'
+    pp_coefs = '1.0'
+    b = ${fparse -p_initial}
+  []
+[]
+
+[Outputs]
+  exodus = true
+  print_linear_residuals = false
+[]
diff --git a/modules/navier_stokes/test/tests/finite_volume/materials/ergun/ergun_in.e b/modules/navier_stokes/test/tests/finite_volume/materials/ergun/ergun_in.e
new file mode 100644
index 000000000000..8417cded883e
Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/materials/ergun/ergun_in.e differ
diff --git a/modules/navier_stokes/test/tests/finite_volume/materials/ergun/gold/ergun_out.e b/modules/navier_stokes/test/tests/finite_volume/materials/ergun/gold/ergun_out.e
new file mode 100644
index 000000000000..4a4aca20ce4f
Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/materials/ergun/gold/ergun_out.e differ
diff --git a/modules/navier_stokes/test/tests/finite_volume/materials/ergun/tests b/modules/navier_stokes/test/tests/finite_volume/materials/ergun/tests
new file mode 100644
index 000000000000..6923332c1497
--- /dev/null
+++ b/modules/navier_stokes/test/tests/finite_volume/materials/ergun/tests
@@ -0,0 +1,10 @@
+[Tests]
+  [ergun]
+    type = Exodiff
+    input = 'ergun.i'
+    exodiff = 'ergun_out.e'
+    requirement = "The system shall be able to model flow around a bend with the porous incompressible Navier Stokes equations using a finite volume discretization and an Ergun drag correlation."
+    issues = '#16756'
+    design = 'FunctorErgunDragCoefficients.md'
+  []
+[]
diff --git a/modules/navier_stokes/test/tests/finite_volume/materials/flow_diode/friction.i b/modules/navier_stokes/test/tests/finite_volume/materials/flow_diode/friction.i
index a27aac629378..45f7da5c6edb 100644
--- a/modules/navier_stokes/test/tests/finite_volume/materials/flow_diode/friction.i
+++ b/modules/navier_stokes/test/tests/finite_volume/materials/flow_diode/friction.i
@@ -79,13 +79,13 @@ rho = 10
   [base_friction]
     type = ADGenericVectorFunctorMaterial
     prop_names = 'Darcy Forchheimer'
-    prop_values = '1 1 1 0.1 0.2 0.3'
+    prop_values = '220 240 260 0 0 0'
   []
   [diode]
     type = NSFVFrictionFlowDiodeMaterial
     direction = '1 0 0'
-    additional_linear_resistance = '10 0 0' #'100 10 20'
-    additional_quadratic_resistance = '10 0 0' #'10 4 4'
+    additional_linear_resistance = '4000 0 0'
+    additional_quadratic_resistance = '0 0 0'
     base_linear_friction_coefs = 'Darcy'
     base_quadratic_friction_coefs = 'Forchheimer'
     sum_linear_friction_name = 'diode_linear'
diff --git a/modules/navier_stokes/test/tests/finite_volume/materials/flow_diode/transient_operation.i b/modules/navier_stokes/test/tests/finite_volume/materials/flow_diode/transient_operation.i
index 5507eed7eaad..649bdccbc15e 100644
--- a/modules/navier_stokes/test/tests/finite_volume/materials/flow_diode/transient_operation.i
+++ b/modules/navier_stokes/test/tests/finite_volume/materials/flow_diode/transient_operation.i
@@ -145,7 +145,7 @@ ny = 5
   [base_friction]
     type = ADGenericVectorFunctorMaterial
     prop_names = 'Darcy Forchheimer'
-    prop_values = '1 1 1 0.1 0.2 0.3'
+    prop_values = '220 240 260 0 0 0'
   []
 
   # Material definitions needed for the diode
@@ -153,7 +153,7 @@ ny = 5
     type = NSFVFrictionFlowDiodeMaterial
     # Friction only in X direction
     direction = '-1 0 0'
-    additional_linear_resistance = '100 0 0'
+    additional_linear_resistance = '20000 0 0'
     additional_quadratic_resistance = '0 0 0'
     base_linear_friction_coefs = 'Darcy'
     base_quadratic_friction_coefs = 'Forchheimer'
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/2d-rc-friction.i b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/2d-rc-friction.i
index 47fb23800e11..3b5961d2163d 100644
--- a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/2d-rc-friction.i
+++ b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/2d-rc-friction.i
@@ -94,10 +94,11 @@ velocity_interp_method = 'rc'
     type = PINSFVMomentumFriction
     variable = superficial_vel_x
     momentum_component = 'x'
-    porosity = porosity
     Darcy_name = 'Darcy_coefficient'
     Forchheimer_name = 'Forchheimer_coefficient'
+    mu = ${mu}
     rho = ${rho}
+    speed = speed
   []
 
   [v_advection]
@@ -127,10 +128,11 @@ velocity_interp_method = 'rc'
     type = PINSFVMomentumFriction
     variable = superficial_vel_y
     momentum_component = 'y'
-    porosity = porosity
     Darcy_name = 'Darcy_coefficient'
     Forchheimer_name = 'Forchheimer_coefficient'
     rho = ${rho}
+    speed = speed
+    mu = ${mu}
   []
 
   [mean-pressure]
@@ -212,12 +214,18 @@ velocity_interp_method = 'rc'
   []
 []
 
-[Materials]
+[FunctorMaterials]
   [darcy]
     type = ADGenericVectorFunctorMaterial
     prop_names = 'Darcy_coefficient Forchheimer_coefficient'
     prop_values = '0.1 0.1 0.1 0.1 0.1 0.1'
   []
+  [speec]
+    type = PINSFVSpeedFunctorMaterial
+    superficial_vel_x = superficial_vel_x
+    superficial_vel_y = superficial_vel_y
+    porosity = porosity
+  []
 []
 
 [Executioner]
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/2d-rc-rz-by-parts.i b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/2d-rc-rz-by-parts.i
index 188a07502f83..19ad9dfaaeeb 100644
--- a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/2d-rc-rz-by-parts.i
+++ b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/2d-rc-rz-by-parts.i
@@ -14,9 +14,6 @@ velocity_interp_method = 'rc'
     nx = 40
     ny = 10
   []
-[]
-
-[Problem]
   coord_type = 'RZ'
   rz_coord_axis = 'X'
 []
@@ -95,10 +92,11 @@ velocity_interp_method = 'rc'
     type = PINSFVMomentumFriction
     variable = u
     momentum_component = 'x'
-    porosity = porosity
     Darcy_name = 'Darcy_coefficient'
     Forchheimer_name = 'Forchheimer_coefficient'
     rho = ${rho}
+    speed = speed
+    mu = ${mu}
   []
 
   [v_advection]
@@ -142,10 +140,11 @@ velocity_interp_method = 'rc'
     type = PINSFVMomentumFriction
     variable = v
     momentum_component = 'y'
-    porosity = porosity
     Darcy_name = 'Darcy_coefficient'
     Forchheimer_name = 'Forchheimer_coefficient'
     rho = ${rho}
+    speed = speed
+    mu = ${mu}
   []
 []
 
@@ -220,12 +219,18 @@ velocity_interp_method = 'rc'
   []
 []
 
-[Materials]
+[FunctorMaterials]
   [darcy]
     type = ADGenericVectorFunctorMaterial
     prop_names = 'Darcy_coefficient Forchheimer_coefficient'
     prop_values = '0.1 0.1 0.1 0.1 0.1 0.1'
   []
+  [speed]
+    type = PINSFVSpeedFunctorMaterial
+    superficial_vel_x = u
+    superficial_vel_y = v
+    porosity = porosity
+  []
 []
 
 [Executioner]
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/2d-rc-rz-by-parts_out.e b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/2d-rc-rz-by-parts_out.e
index d4c3aa2923e2..235c4e6d20d6 100644
Binary files a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/2d-rc-rz-by-parts_out.e and b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/2d-rc-rz-by-parts_out.e differ
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/linear_friction_out.e b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/linear_friction_out.e
index fa17d7b54839..710c89a31298 100644
Binary files a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/linear_friction_out.e and b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/linear_friction_out.e differ
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction-block.e b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction-block.e
index c345b00a8a86..66c191523b80 100644
Binary files a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction-block.e and b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction-block.e differ
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction-free-slip.e b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction-free-slip.e
index 804958f62c45..bd75faa16955 100644
Binary files a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction-free-slip.e and b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction-free-slip.e differ
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction.e b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction.e
index dd977d40a421..39702ef62df6 100644
Binary files a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction.e and b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/gold/rc-friction.e differ
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/segregated/2d-momentum.i b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/segregated/2d-momentum.i
index cdf923f556bf..f0357051a13a 100644
--- a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/segregated/2d-momentum.i
+++ b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/segregated/2d-momentum.i
@@ -94,8 +94,9 @@ pressure_tag = "pressure_grad"
     momentum_component = 'y'
     Darcy_name = 'Darcy_coefficient'
     Forchheimer_name = 'Forchheimer_coefficient'
-    porosity = 'porosity'
     rho = ${rho}
+    speed = speed
+    mu = ${mu}
   []
 
   [v_advection]
@@ -126,8 +127,9 @@ pressure_tag = "pressure_grad"
     momentum_component = 'y'
     Darcy_name = 'Darcy_coefficient'
     Forchheimer_name = 'Forchheimer_coefficient'
-    porosity = 'porosity'
     rho = ${rho}
+    speed = speed
+    mu = ${mu}
   []
 
   [p_diffusion]
@@ -215,12 +217,18 @@ pressure_tag = "pressure_grad"
   #####################################################################
 []
 
-[Materials]
+[FunctorMaterials]
   [darcy]
     type = ADGenericVectorFunctorMaterial
     prop_names = 'Darcy_coefficient Forchheimer_coefficient'
     prop_values = '0.01 0.02 0.03 0.01 0.02 0.03'
   []
+  [speed]
+    type = PINSFVSpeedFunctorMaterial
+    superficial_vel_x = superficial_vel_x
+    superficial_vel_y = superficial_vel_y
+    porosity = porosity
+  []
 []
 
 [Executioner]
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/segregated/gold/2d-momentum-slip-friction.e b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/segregated/gold/2d-momentum-slip-friction.e
index 4eafe0ce2eec..a5a9ef23251e 100644
Binary files a/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/segregated/gold/2d-momentum-slip-friction.e and b/modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/segregated/gold/2d-momentum-slip-friction.e differ
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-corrected.py b/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-corrected.py
index 1dafe93cec2e..5b0d6a8801c3 100755
--- a/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-corrected.py
+++ b/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-corrected.py
@@ -11,8 +11,8 @@
 
 porosity = '.5 + .1 * sin(pi * x / 4) * cos(pi * y / 4)'
 
-f_u, e_u = mms.evaluate('div(vel*rho*u/porosity) - div(mu*porosity*grad(u/porosity)) + porosity*grad(p).dot(e_i) + (darcy + forch)*rho*u/porosity', u, variable='u', vel=vel, p=p, porosity=porosity, scalars=['mu', 'rho', 'darcy', 'forch'])
-f_v, e_v = mms.evaluate('div(vel*rho*v/porosity) - div(mu*porosity*grad(v/porosity)) + porosity*grad(p).dot(e_j) + (darcy + forch)*rho*v/porosity', v, variable='v', vel=vel, p=p, porosity=porosity, scalars=['mu', 'rho', 'darcy', 'forch'])
+f_u, e_u = mms.evaluate('div(vel*rho*u/porosity) - div(mu*porosity*grad(u/porosity)) + porosity*grad(p).dot(e_i) + darcy*mu*u + forch*rho/2*u*sqrt(u*u+v*v)/porosity', u, variable='u', v=v, vel=vel, p=p, porosity=porosity, scalars=['mu', 'rho', 'darcy', 'forch'])
+f_v, e_v = mms.evaluate('div(vel*rho*v/porosity) - div(mu*porosity*grad(v/porosity)) + porosity*grad(p).dot(e_j) + darcy*mu*v + forch*rho/2*v*sqrt(u*u+v*v)/porosity', v, variable='v', u=u, vel=vel, p=p, porosity=porosity, scalars=['mu', 'rho', 'darcy', 'forch'])
 f_p, e_p = mms.evaluate('div(vel*rho)', p, variable='p', vel=vel, scalars=['rho'])
 
 rho = sympy.Symbol('rho')
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-interpolation-corrected-action.i b/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-interpolation-corrected-action.i
index 14bbe67d0837..d9ef505a3585 100644
--- a/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-interpolation-corrected-action.i
+++ b/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-interpolation-corrected-action.i
@@ -110,7 +110,7 @@ forch=1.1
   []
   [forcing_u]
     type = ParsedFunction
-    expression = '-mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(-1/4*pi^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.025*pi^2*sin((1/4)*x*pi)*sin((1/4)*y*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00625*pi^2*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00125*pi^2*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^3) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(-1/4*pi^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.00625*pi^2*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.025*pi^2*sin((1/2)*x*pi)*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/4)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00125*pi^2*sin((1/2)*y*pi)*cos((1/4)*x*pi)^2*cos((1/2)*x*pi)*cos((1/4)*y*pi)^2/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^3) - 0.025*pi*mu*(-1/2*pi*sin((1/2)*x*pi)*sin((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 0.025*pi*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2)*cos((1/4)*x*pi)*cos((1/4)*y*pi) + 0.025*pi*mu*((1/2)*pi*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.025*pi*sin((1/4)*x*pi)*sin((1/4)*y*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2)*sin((1/4)*x*pi)*sin((1/4)*y*pi) + rho*(darcy + forch)*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 1/2*pi*rho*sin((1/4)*x*pi)*sin((1/2)*y*pi)^2*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + (1/2)*pi*rho*sin((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)^2/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - pi*rho*sin((1/2)*x*pi)*sin((1/2)*y*pi)^2*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.025*pi*rho*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 - 0.025*pi*rho*sin((1/2)*y*pi)^2*cos((1/4)*x*pi)*cos((1/2)*x*pi)^2*cos((1/4)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 - 1/4*pi*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*sin((1/4)*x*pi)*sin((3/2)*y*pi)'
+    expression = 'darcy*mu*sin((1/2)*y*pi)*cos((1/2)*x*pi) + (1/2)*forch*rho*sqrt(sin((1/4)*x*pi)^2*cos((1/2)*y*pi)^2 + sin((1/2)*y*pi)^2*cos((1/2)*x*pi)^2)*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(0.1*pi^2*sin((1/4)*x*pi)*sin((1/4)*y*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.025*pi^2*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.01*pi^2*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^3 - 1/4*pi^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(0.025*pi^2*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.1*pi^2*sin((1/2)*x*pi)*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/4)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.01*pi^2*sin((1/2)*y*pi)*cos((1/4)*x*pi)^2*cos((1/2)*x*pi)*cos((1/4)*y*pi)^2/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^3 - 1/4*pi^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)) + 0.025*pi*mu*(0.1*pi*sin((1/4)*x*pi)*sin((1/4)*y*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + (1/2)*pi*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5))*sin((1/4)*x*pi)*sin((1/4)*y*pi) - 0.025*pi*mu*(-0.1*pi*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 1/2*pi*sin((1/2)*x*pi)*sin((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5))*cos((1/4)*x*pi)*cos((1/4)*y*pi) + 0.1*pi*rho*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 0.1*pi*rho*sin((1/2)*y*pi)^2*cos((1/4)*x*pi)*cos((1/2)*x*pi)^2*cos((1/4)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 1/2*pi*rho*sin((1/4)*x*pi)*sin((1/2)*y*pi)^2*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + (1/2)*pi*rho*sin((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)^2/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - pi*rho*sin((1/2)*x*pi)*sin((1/2)*y*pi)^2*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 1/4*pi*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*sin((1/4)*x*pi)*sin((3/2)*y*pi)'
     symbol_names = 'mu rho darcy forch'
     symbol_values = '${mu} ${rho} ${darcy} ${forch}'
   []
@@ -120,7 +120,7 @@ forch=1.1
   []
   [forcing_v]
     type = ParsedFunction
-    expression = '-mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(-1/4*pi^2*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 0.025*pi^2*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*sin((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00625*pi^2*sin((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00125*pi^2*sin((1/4)*x*pi)^3*sin((1/4)*y*pi)^2*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^3) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(-1/16*pi^2*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.00625*pi^2*sin((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 - 0.0125*pi^2*cos((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00125*pi^2*sin((1/4)*x*pi)*cos((1/4)*x*pi)^2*cos((1/4)*y*pi)^2*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^3) + 0.025*pi*mu*(-1/2*pi*sin((1/4)*x*pi)*sin((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.025*pi*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2)*sin((1/4)*x*pi)*sin((1/4)*y*pi) - 0.025*pi*mu*((1/4)*pi*cos((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 0.025*pi*sin((1/4)*x*pi)*cos((1/4)*x*pi)*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2)*cos((1/4)*x*pi)*cos((1/4)*y*pi) + rho*(darcy + forch)*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - pi*rho*sin((1/4)*x*pi)^2*sin((1/2)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 1/2*pi*rho*sin((1/4)*x*pi)*sin((1/2)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + (1/4)*pi*rho*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.025*pi*rho*sin((1/4)*x*pi)^3*sin((1/4)*y*pi)*cos((1/2)*y*pi)^2/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 - 0.025*pi*rho*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + (3/2)*pi*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*cos((1/4)*x*pi)*cos((3/2)*y*pi)'
+    expression = 'darcy*mu*sin((1/4)*x*pi)*cos((1/2)*y*pi) + (1/2)*forch*rho*sqrt(sin((1/4)*x*pi)^2*cos((1/2)*y*pi)^2 + sin((1/2)*y*pi)^2*cos((1/2)*x*pi)^2)*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(-0.1*pi^2*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*sin((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.025*pi^2*sin((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.01*pi^2*sin((1/4)*x*pi)^3*sin((1/4)*y*pi)^2*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^3 - 1/4*pi^2*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(0.025*pi^2*sin((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 0.05*pi^2*cos((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.01*pi^2*sin((1/4)*x*pi)*cos((1/4)*x*pi)^2*cos((1/4)*y*pi)^2*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^3 - 1/16*pi^2*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)) + 0.025*pi*mu*(0.1*pi*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 1/2*pi*sin((1/4)*x*pi)*sin((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5))*sin((1/4)*x*pi)*sin((1/4)*y*pi) - 0.025*pi*mu*(-0.1*pi*sin((1/4)*x*pi)*cos((1/4)*x*pi)*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + (1/4)*pi*cos((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5))*cos((1/4)*x*pi)*cos((1/4)*y*pi) + 0.1*pi*rho*sin((1/4)*x*pi)^3*sin((1/4)*y*pi)*cos((1/2)*y*pi)^2/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 0.1*pi*rho*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - pi*rho*sin((1/4)*x*pi)^2*sin((1/2)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 1/2*pi*rho*sin((1/4)*x*pi)*sin((1/2)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + (1/4)*pi*rho*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + (3/2)*pi*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*cos((1/4)*x*pi)*cos((3/2)*y*pi)'
     symbol_names = 'mu rho darcy forch'
     symbol_values = '${mu} ${rho} ${darcy} ${forch}'
   []
diff --git a/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-interpolation-corrected.i b/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-interpolation-corrected.i
index 4f0cbb62f3aa..5dea67676d19 100644
--- a/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-interpolation-corrected.i
+++ b/modules/navier_stokes/test/tests/finite_volume/pins/mms/porosity_change/pressure-interpolation-corrected.i
@@ -105,15 +105,17 @@ forch = 1.1
     type = PINSFVMomentumFriction
     variable = u
     momentum_component = 'x'
-    porosity = porosity
     rho = ${rho}
+    speed = speed
+    mu = ${mu}
   []
   [u_correction]
     type = PINSFVMomentumFrictionCorrection
     variable = u
     momentum_component = 'x'
-    porosity = porosity
     rho = ${rho}
+    speed = speed
+    mu = ${mu}
   []
   [u_forcing]
     type = INSFVBodyForce
@@ -148,15 +150,17 @@ forch = 1.1
     type = PINSFVMomentumFriction
     variable = v
     momentum_component = 'y'
-    porosity = porosity
     rho = ${rho}
+    speed = speed
+    mu = ${mu}
   []
   [v_correction]
     type = PINSFVMomentumFrictionCorrection
     variable = v
     momentum_component = 'y'
-    porosity = porosity
     rho = ${rho}
+    speed = speed
+    mu = ${mu}
   []
   [v_forcing]
     type = INSFVBodyForce
@@ -199,12 +203,18 @@ forch = 1.1
   []
 []
 
-[Materials]
+[FunctorMaterials]
   [darcy]
     type = ADGenericVectorFunctorMaterial
     prop_names = 'Darcy_coefficient Forchheimer_coefficient'
     prop_values = '${darcy} ${darcy} ${darcy} ${forch} ${forch} ${forch}'
   []
+  [speed]
+    type = PINSFVSpeedFunctorMaterial
+    superficial_vel_x = u
+    superficial_vel_y = v
+    porosity = porosity
+  []
 []
 
 [Functions]
@@ -219,7 +229,7 @@ forch = 1.1
   []
   [forcing_u]
     type = ParsedFunction
-    expression = '-mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(-1/4*pi^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.025*pi^2*sin((1/4)*x*pi)*sin((1/4)*y*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00625*pi^2*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00125*pi^2*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^3) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(-1/4*pi^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.00625*pi^2*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.025*pi^2*sin((1/2)*x*pi)*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/4)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00125*pi^2*sin((1/2)*y*pi)*cos((1/4)*x*pi)^2*cos((1/2)*x*pi)*cos((1/4)*y*pi)^2/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^3) - 0.025*pi*mu*(-1/2*pi*sin((1/2)*x*pi)*sin((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 0.025*pi*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2)*cos((1/4)*x*pi)*cos((1/4)*y*pi) + 0.025*pi*mu*((1/2)*pi*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.025*pi*sin((1/4)*x*pi)*sin((1/4)*y*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2)*sin((1/4)*x*pi)*sin((1/4)*y*pi) + rho*(darcy + forch)*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 1/2*pi*rho*sin((1/4)*x*pi)*sin((1/2)*y*pi)^2*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + (1/2)*pi*rho*sin((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)^2/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - pi*rho*sin((1/2)*x*pi)*sin((1/2)*y*pi)^2*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.025*pi*rho*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 - 0.025*pi*rho*sin((1/2)*y*pi)^2*cos((1/4)*x*pi)*cos((1/2)*x*pi)^2*cos((1/4)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 - 1/4*pi*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*sin((1/4)*x*pi)*sin((3/2)*y*pi)'
+    expression = 'darcy*mu*sin((1/2)*y*pi)*cos((1/2)*x*pi) + (1/2)*forch*rho*sqrt(sin((1/4)*x*pi)^2*cos((1/2)*y*pi)^2 + sin((1/2)*y*pi)^2*cos((1/2)*x*pi)^2)*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(0.1*pi^2*sin((1/4)*x*pi)*sin((1/4)*y*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.025*pi^2*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.01*pi^2*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^3 - 1/4*pi^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(0.025*pi^2*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.1*pi^2*sin((1/2)*x*pi)*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/4)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.01*pi^2*sin((1/2)*y*pi)*cos((1/4)*x*pi)^2*cos((1/2)*x*pi)*cos((1/4)*y*pi)^2/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^3 - 1/4*pi^2*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)) + 0.025*pi*mu*(0.1*pi*sin((1/4)*x*pi)*sin((1/4)*y*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + (1/2)*pi*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5))*sin((1/4)*x*pi)*sin((1/4)*y*pi) - 0.025*pi*mu*(-0.1*pi*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 1/2*pi*sin((1/2)*x*pi)*sin((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5))*cos((1/4)*x*pi)*cos((1/4)*y*pi) + 0.1*pi*rho*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*sin((1/2)*y*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 0.1*pi*rho*sin((1/2)*y*pi)^2*cos((1/4)*x*pi)*cos((1/2)*x*pi)^2*cos((1/4)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 1/2*pi*rho*sin((1/4)*x*pi)*sin((1/2)*y*pi)^2*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + (1/2)*pi*rho*sin((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)^2/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - pi*rho*sin((1/2)*x*pi)*sin((1/2)*y*pi)^2*cos((1/2)*x*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 1/4*pi*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*sin((1/4)*x*pi)*sin((3/2)*y*pi)'
     symbol_names = 'mu rho darcy forch'
     symbol_values = '${mu} ${rho} ${darcy} ${forch}'
   []
@@ -229,7 +239,7 @@ forch = 1.1
   []
   [forcing_v]
     type = ParsedFunction
-    expression = '-mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(-1/4*pi^2*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 0.025*pi^2*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*sin((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00625*pi^2*sin((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00125*pi^2*sin((1/4)*x*pi)^3*sin((1/4)*y*pi)^2*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^3) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(-1/16*pi^2*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.00625*pi^2*sin((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 - 0.0125*pi^2*cos((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + 0.00125*pi^2*sin((1/4)*x*pi)*cos((1/4)*x*pi)^2*cos((1/4)*y*pi)^2*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^3) + 0.025*pi*mu*(-1/2*pi*sin((1/4)*x*pi)*sin((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.025*pi*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2)*sin((1/4)*x*pi)*sin((1/4)*y*pi) - 0.025*pi*mu*((1/4)*pi*cos((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 0.025*pi*sin((1/4)*x*pi)*cos((1/4)*x*pi)*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2)*cos((1/4)*x*pi)*cos((1/4)*y*pi) + rho*(darcy + forch)*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - pi*rho*sin((1/4)*x*pi)^2*sin((1/2)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 1/2*pi*rho*sin((1/4)*x*pi)*sin((1/2)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + (1/4)*pi*rho*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + 0.025*pi*rho*sin((1/4)*x*pi)^3*sin((1/4)*y*pi)*cos((1/2)*y*pi)^2/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 - 0.025*pi*rho*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)^2 + (3/2)*pi*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*cos((1/4)*x*pi)*cos((3/2)*y*pi)'
+    expression = 'darcy*mu*sin((1/4)*x*pi)*cos((1/2)*y*pi) + (1/2)*forch*rho*sqrt(sin((1/4)*x*pi)^2*cos((1/2)*y*pi)^2 + sin((1/2)*y*pi)^2*cos((1/2)*x*pi)^2)*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(-0.1*pi^2*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*sin((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.025*pi^2*sin((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.01*pi^2*sin((1/4)*x*pi)^3*sin((1/4)*y*pi)^2*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^3 - 1/4*pi^2*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)) - mu*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*(0.025*pi^2*sin((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 0.05*pi^2*cos((1/4)*x*pi)^2*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + 0.01*pi^2*sin((1/4)*x*pi)*cos((1/4)*x*pi)^2*cos((1/4)*y*pi)^2*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^3 - 1/16*pi^2*sin((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)) + 0.025*pi*mu*(0.1*pi*sin((1/4)*x*pi)^2*sin((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 1/2*pi*sin((1/4)*x*pi)*sin((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5))*sin((1/4)*x*pi)*sin((1/4)*y*pi) - 0.025*pi*mu*(-0.1*pi*sin((1/4)*x*pi)*cos((1/4)*x*pi)*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 + (1/4)*pi*cos((1/4)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5))*cos((1/4)*x*pi)*cos((1/4)*y*pi) + 0.1*pi*rho*sin((1/4)*x*pi)^3*sin((1/4)*y*pi)*cos((1/2)*y*pi)^2/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - 0.1*pi*rho*sin((1/4)*x*pi)*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/4)*y*pi)*cos((1/2)*y*pi)/(0.2*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 1)^2 - pi*rho*sin((1/4)*x*pi)^2*sin((1/2)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) - 1/2*pi*rho*sin((1/4)*x*pi)*sin((1/2)*x*pi)*sin((1/2)*y*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + (1/4)*pi*rho*sin((1/2)*y*pi)*cos((1/4)*x*pi)*cos((1/2)*x*pi)*cos((1/2)*y*pi)/(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5) + (3/2)*pi*(0.1*sin((1/4)*x*pi)*cos((1/4)*y*pi) + 0.5)*cos((1/4)*x*pi)*cos((3/2)*y*pi)'
     symbol_names = 'mu rho darcy forch'
     symbol_values = '${mu} ${rho} ${darcy} ${forch}'
   []
diff --git a/modules/navier_stokes/test/tests/finite_volume/pwcns/channel-flow/2d-transient-gas.i b/modules/navier_stokes/test/tests/finite_volume/pwcns/channel-flow/2d-transient-gas.i
index 4a67669bcd19..3e7bb7f25dec 100644
--- a/modules/navier_stokes/test/tests/finite_volume/pwcns/channel-flow/2d-transient-gas.i
+++ b/modules/navier_stokes/test/tests/finite_volume/pwcns/channel-flow/2d-transient-gas.i
@@ -312,7 +312,6 @@ velocity_interp_method = 'rc'
   [speed]
     type = PINSFVSpeedFunctorMaterial
     porosity = 'porosity'
-    T_fluid = 'T_fluid'
     superficial_vel_x = 'superficial_vel_x'
     superficial_vel_y = 'superficial_vel_y'
   []
diff --git a/modules/navier_stokes/test/tests/finite_volume/wcns/natural_convection/natural_circulation_pipe.i b/modules/navier_stokes/test/tests/finite_volume/wcns/natural_convection/natural_circulation_pipe.i
index e73ec04ad979..f88414fea4ac 100644
--- a/modules/navier_stokes/test/tests/finite_volume/wcns/natural_convection/natural_circulation_pipe.i
+++ b/modules/navier_stokes/test/tests/finite_volume/wcns/natural_convection/natural_circulation_pipe.i
@@ -73,14 +73,27 @@ gamma = 1.4
     pressure_function = '${fparse p0}'
     momentum_advection_interpolation = 'upwind'
     mass_advection_interpolation = 'upwind'
-    friction_types = 'Darcy'
-    friction_coeffs = 'Darcy_coef'
     porous_medium_treatment = true
     porosity = porosity
     energy_advection_interpolation = 'average'
   []
 []
 
+[FVKernels]
+  [u_friction]
+    type = INSFVMomentumFriction
+    variable = superficial_vel_x
+    linear_coef_name = linear_friction_coeff
+    momentum_component = 'x'
+  []
+  [v_friction]
+    type = INSFVMomentumFriction
+    variable = superficial_vel_y
+    linear_coef_name = linear_friction_coeff
+    momentum_component = 'y'
+  []
+[]
+
 [Executioner]
   type = Transient
   solve_type = 'NEWTON'
@@ -106,7 +119,7 @@ gamma = 1.4
   []
 []
 
-[Materials]
+[FunctorMaterials]
   [fluid_props_to_mat_props]
     type = GeneralFunctorFluidProps
     fp = air
@@ -121,14 +134,15 @@ gamma = 1.4
 
   [scalar_props]
     type = ADGenericFunctorMaterial
-    prop_names = 'porosity'
-    prop_values = '1      '
+    prop_names = 'porosity loss_coeff'
+    prop_values = '1       1.3'
   []
 
-  [vector_props]
-    type = ADGenericVectorFunctorMaterial
-    prop_names = 'Darcy_coef'
-    prop_values = '1.3 1.3 1.3'
+  [linear_friction_coeff]
+    type = ADParsedFunctorMaterial
+    property_name = 'linear_friction_coeff'
+    expression = 'loss_coeff * rho'
+    functor_names = 'loss_coeff rho'
   []
 []