diff --git a/.coverage b/.coverage
index 2f6f19017d77..67a7614684dc 100644
--- a/.coverage
+++ b/.coverage
@@ -36,7 +36,7 @@ require_total = 85
 require_total = 84
 
 [fsi]
-require_total = 90
+require_total = 85
 
 [functional_expansion_tools]
 require_total = 82
diff --git a/framework/include/interfacekernels/ADInterfaceKernel.h b/framework/include/interfacekernels/ADInterfaceKernel.h
index bfcce5c3ba12..208fe21c4cce 100644
--- a/framework/include/interfacekernels/ADInterfaceKernel.h
+++ b/framework/include/interfacekernels/ADInterfaceKernel.h
@@ -35,28 +35,27 @@ class ADInterfaceKernelTempl : public InterfaceKernelBase, public NeighborMooseV
    * Using the passed DGResidual type, selects the correct test function space and residual block,
    * and then calls computeQpResidual
    */
-  void computeElemNeighResidual(Moose::DGResidualType type) override final;
+  void computeElemNeighResidual(Moose::DGResidualType type);
 
   /**
    * Using the passed DGJacobian type, selects the correct test function and trial function spaces
    * and
    * jacobian block, and then calls computeQpJacobian
    */
-  void computeElemNeighJacobian(Moose::DGJacobianType type) override final;
+  void computeElemNeighJacobian(Moose::DGJacobianType type);
 
   /**
    * Using the passed DGJacobian type, selects the correct test function and trial function spaces
    * and
    * jacobian block, and then calls computeQpOffDiagJacobian with the passed jvar
    */
-  void computeOffDiagElemNeighJacobian(Moose::DGJacobianType type,
-                                       unsigned int jvar) override final;
+  void computeOffDiagElemNeighJacobian(Moose::DGJacobianType type, unsigned int jvar);
 
   /// Selects the correct Jacobian type and routine to call for the primary variable jacobian
-  void computeElementOffDiagJacobian(unsigned int jvar) override final;
+  virtual void computeElementOffDiagJacobian(unsigned int jvar) override final;
 
   /// Selects the correct Jacobian type and routine to call for the secondary variable jacobian
-  void computeNeighborOffDiagJacobian(unsigned int jvar) override final;
+  virtual void computeNeighborOffDiagJacobian(unsigned int jvar) override final;
 
   /// Computes the residual for the current side.
   void computeResidual() override final;
diff --git a/framework/include/interfacekernels/InterfaceKernel.h b/framework/include/interfacekernels/InterfaceKernel.h
index 5093628c5d2a..071186ef23c9 100644
--- a/framework/include/interfacekernels/InterfaceKernel.h
+++ b/framework/include/interfacekernels/InterfaceKernel.h
@@ -48,22 +48,21 @@ class InterfaceKernelTempl : public InterfaceKernelBase, public NeighborMooseVar
    * Using the passed DGResidual type, selects the correct test function space and residual block,
    * and then calls computeQpResidual
    */
-  virtual void computeElemNeighResidual(Moose::DGResidualType type) override;
+  virtual void computeElemNeighResidual(Moose::DGResidualType type);
 
   /**
    * Using the passed DGJacobian type, selects the correct test function and trial function spaces
    * and
    * jacobian block, and then calls computeQpJacobian
    */
-  virtual void computeElemNeighJacobian(Moose::DGJacobianType type) override;
+  virtual void computeElemNeighJacobian(Moose::DGJacobianType type);
 
   /**
    * Using the passed DGJacobian type, selects the correct test function and trial function spaces
    * and
    * jacobian block, and then calls computeQpOffDiagJacobian with the passed jvar
    */
-  virtual void computeOffDiagElemNeighJacobian(Moose::DGJacobianType type,
-                                               unsigned int jvar) override;
+  virtual void computeOffDiagElemNeighJacobian(Moose::DGJacobianType type, unsigned int jvar);
 
   /// Selects the correct Jacobian type and routine to call for the primary variable jacobian
   virtual void computeElementOffDiagJacobian(unsigned int jvar) override;
diff --git a/framework/include/interfacekernels/InterfaceKernelBase.h b/framework/include/interfacekernels/InterfaceKernelBase.h
index af70367dd128..c499cac15f2b 100644
--- a/framework/include/interfacekernels/InterfaceKernelBase.h
+++ b/framework/include/interfacekernels/InterfaceKernelBase.h
@@ -35,26 +35,6 @@ class InterfaceKernelBase : public NeighborResidualObject,
   /// The neighbor variable number that this interface kernel operates on
   virtual const MooseVariableFEBase & neighborVariable() const = 0;
 
-  /**
-   * Using the passed DGResidual type, selects the correct test function space and residual block,
-   * and then calls computeQpResidual
-   */
-  virtual void computeElemNeighResidual(Moose::DGResidualType type) = 0;
-
-  /**
-   * Using the passed DGJacobian type, selects the correct test function and trial function spaces
-   * and
-   * jacobian block, and then calls computeQpJacobian
-   */
-  virtual void computeElemNeighJacobian(Moose::DGJacobianType type) = 0;
-
-  /**
-   * Using the passed DGJacobian type, selects the correct test function and trial function spaces
-   * and
-   * jacobian block, and then calls computeQpOffDiagJacobian with the passed jvar
-   */
-  virtual void computeOffDiagElemNeighJacobian(Moose::DGJacobianType type, unsigned int jvar) = 0;
-
   /// Selects the correct Jacobian type and routine to call for the primary variable jacobian
   virtual void computeElementOffDiagJacobian(unsigned int jvar) = 0;
 
diff --git a/framework/src/kernels/ADKernelSUPG.C b/framework/src/kernels/ADKernelSUPG.C
index 00f41a1c7c09..ec397112f765 100644
--- a/framework/src/kernels/ADKernelSUPG.C
+++ b/framework/src/kernels/ADKernelSUPG.C
@@ -21,7 +21,9 @@ ADKernelSUPGTempl<T>::validParams()
   InputParameters params = ADKernelStabilizedTempl<T>::validParams();
   params.addParam<MaterialPropertyName>(
       "tau_name", "tau", "The name of the stabilization parameter tau.");
-  params.addRequiredCoupledVar("velocity", "The velocity variable.");
+  params.addCoupledVar("velocity", "The velocity variable.");
+  params.addParam<MaterialPropertyName>("material_velocity",
+                                        "A material property describing the velocity");
   return params;
 }
 
@@ -29,7 +31,10 @@ template <typename T>
 ADKernelSUPGTempl<T>::ADKernelSUPGTempl(const InputParameters & parameters)
   : ADKernelStabilizedTempl<T>(parameters),
     _tau(this->template getADMaterialProperty<Real>("tau_name")),
-    _velocity(this->adCoupledVectorValue("velocity"))
+    _velocity(
+        this->isParamValid("velocity")
+            ? this->adCoupledVectorValue("velocity")
+            : this->template getADMaterialProperty<RealVectorValue>("material_velocity").get())
 {
 }
 
diff --git a/modules/fsi/doc/content/source/interfacekernels/ADPenaltyVelocityContinuity.md b/modules/fsi/doc/content/source/interfacekernels/ADPenaltyVelocityContinuity.md
new file mode 100644
index 000000000000..8c493d9c47ea
--- /dev/null
+++ b/modules/fsi/doc/content/source/interfacekernels/ADPenaltyVelocityContinuity.md
@@ -0,0 +1,15 @@
+# ADPenaltyVelocityContinuity
+
+!syntax description /InterfaceKernels/ADPenaltyVelocityContinuity
+
+This object is meant for use in coupling the Navier-Stokes fluid equations with
+solid mechanics equations at a fluid-structure interface. In that context this
+object imposes both continuity of stress and continuity of material velocity via
+a penalty condition. For more information on the error of penalty methods,
+please see [CoupledPenaltyInterfaceDiffusion.md].
+
+!syntax parameters /InterfaceKernels/ADPenaltyVelocityContinuity
+
+!syntax inputs /InterfaceKernels/ADPenaltyVelocityContinuity
+
+!syntax children /InterfaceKernels/ADPenaltyVelocityContinuity
diff --git a/modules/fsi/doc/content/source/kernels/ConvectedMesh.md b/modules/fsi/doc/content/source/kernels/ConvectedMesh.md
index 6a27801953af..b1218d9ac01a 100644
--- a/modules/fsi/doc/content/source/kernels/ConvectedMesh.md
+++ b/modules/fsi/doc/content/source/kernels/ConvectedMesh.md
@@ -1,11 +1,12 @@
 # ConvectedMesh
 
-## Description
+!syntax description /Kernels/ConvectedMesh
 
 `ConvectedMesh` implements the corresponding weak form for the components of
 the term:
 
 \begin{equation}
+\label{convection}
 -\rho \left(\frac{\partial\vec{d_m}}{\partial t} \cdot \nabla\right) \vec{u}
 \end{equation}
 
@@ -13,3 +14,9 @@ where $\rho$ is the density, $\vec{d_m}$ is the fluid mesh displacements, and
 $\vec{u}$ is the fluid velocity. This term is essential for obtaining the
 correct convective derivative of the fluid in cases where the fluid mesh is
 dynamic, e.g. in simulations of fluid-structure interaction.
+
+!syntax parameters /Kernels/ConvectedMesh
+
+!syntax inputs /Kernels/ConvectedMesh
+
+!syntax children /Kernels/ConvectedMesh
diff --git a/modules/fsi/doc/content/source/kernels/ConvectedMeshPSPG.md b/modules/fsi/doc/content/source/kernels/ConvectedMeshPSPG.md
new file mode 100644
index 000000000000..292f3d99303c
--- /dev/null
+++ b/modules/fsi/doc/content/source/kernels/ConvectedMeshPSPG.md
@@ -0,0 +1,31 @@
+# ConvectedMeshPSPG
+
+!syntax description /Kernels/ConvectedMeshPSPG
+
+This object adds the pressure-stabilized Petrov-Galerkin term to the pressure
+equation corresponding to the [ConvectedMesh.md] object. It implements the
+weak form
+
+\begin{equation}
+\label{pspg}
+\int \nabla \psi \cdot \left(\tau\frac{\partial\vec{d_m}}{\partial t}\cdot\nabla\vec{u}\right) dV
+\end{equation}
+
+where $\nabla \psi$ is the gradient of the pressure test function, $\tau$ is a
+stabilization parameter computed automatically by the `navier_stokes` base class
+`INSBase`, $\vec{d_m}$ is the fluid mesh displacements, and
+$\vec{u}$ is the fluid velocity. Note that when comparing [pspg] with
+[!eqref](ConvectedMesh.md#convection) that the minus sign and density $\rho$ have
+disappeared. This is because the form of PSPG stabilization is
+
+\begin{equation}
+\int \nabla \psi \cdot \left(-\frac{\tau}{\rho}\vec{R}\right)
+\end{equation}
+
+where $\vec{R}$ denotes the strong form of the momentum residuals.
+
+!syntax parameters /Kernels/ConvectedMeshPSPG
+
+!syntax inputs /Kernels/ConvectedMeshPSPG
+
+!syntax children /Kernels/ConvectedMeshPSPG
diff --git a/modules/fsi/include/interfacekernels/ADPenaltyVelocityContinuity.h b/modules/fsi/include/interfacekernels/ADPenaltyVelocityContinuity.h
new file mode 100644
index 000000000000..3611710689f6
--- /dev/null
+++ b/modules/fsi/include/interfacekernels/ADPenaltyVelocityContinuity.h
@@ -0,0 +1,101 @@
+///* 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 "InterfaceKernelBase.h"
+#include "MooseVariableFE.h"
+
+/**
+ * Interface kernel for enforcing continuity of stress and velocity
+ */
+class ADPenaltyVelocityContinuity : public InterfaceKernelBase
+{
+public:
+  static InputParameters validParams();
+
+  ADPenaltyVelocityContinuity(const InputParameters & parameters);
+
+  virtual void computeResidual() override;
+  virtual void computeResidualAndJacobian() override;
+  virtual void computeJacobian() override;
+  virtual void computeElementOffDiagJacobian(unsigned int jvar) override;
+  virtual void computeNeighborOffDiagJacobian(unsigned int jvar) override;
+  virtual const MooseVariableFieldBase & variable() const override;
+  virtual const MooseVariableFieldBase & neighborVariable() const override;
+
+protected:
+  /// The penalty factor
+  const Real _penalty;
+
+  /// Fluid velocity variable
+  const VectorMooseVariable * const _velocity_var;
+
+  /// Fluid velocity values
+  const ADVectorVariableValue & _velocity;
+
+  /// Solid velocity values
+  std::vector<const ADVariableValue *> _solid_velocities;
+
+  /// Displacement variables
+  std::vector<const MooseVariable *> _displacements;
+
+  /// JxW with displacement derivatives
+  const MooseArray<ADReal> & _ad_JxW;
+
+  /// Coordinate transformation with displacement derivatives
+  const MooseArray<ADReal> & _ad_coord;
+
+  /// Residuals data member to avoid constant heap allocation
+  std::vector<ADReal> _residuals;
+
+  /// Jump data member to avoid constant heap allocations
+  std::vector<ADRealVectorValue> _qp_jumps;
+};
+
+inline void
+ADPenaltyVelocityContinuity::computeJacobian()
+{
+  computeResidual();
+}
+
+inline void
+ADPenaltyVelocityContinuity::computeResidualAndJacobian()
+{
+  computeResidual();
+}
+
+inline void
+ADPenaltyVelocityContinuity::computeElementOffDiagJacobian(const unsigned int jvar)
+{
+  if (jvar == _velocity_var->number())
+    // Only need to do this once because AD does everything all at once
+    computeResidual();
+}
+
+inline void
+ADPenaltyVelocityContinuity::computeNeighborOffDiagJacobian(unsigned int)
+{
+}
+
+inline const MooseVariableFieldBase &
+ADPenaltyVelocityContinuity::variable() const
+{
+  return *_velocity_var;
+}
+
+inline const MooseVariableFieldBase &
+ADPenaltyVelocityContinuity::neighborVariable() const
+{
+  if (_displacements.empty() || !_displacements.front())
+    mooseError("The 'neighborVariable' method was called which requires that displacements be "
+               "actual variables.");
+
+  return *_displacements.front();
+}
diff --git a/modules/fsi/include/kernels/ConvectedMesh.h b/modules/fsi/include/kernels/ConvectedMesh.h
index 8b953d79e7de..d288f0ece028 100644
--- a/modules/fsi/include/kernels/ConvectedMesh.h
+++ b/modules/fsi/include/kernels/ConvectedMesh.h
@@ -9,12 +9,12 @@
 
 #pragma once
 
-#include "Kernel.h"
+#include "INSBase.h"
 
 /**
  * This calculates the time derivative for a coupled variable
  **/
-class ConvectedMesh : public Kernel
+class ConvectedMesh : public INSBase
 {
 public:
   static InputParameters validParams();
@@ -26,6 +26,17 @@ class ConvectedMesh : public Kernel
   virtual Real computeQpJacobian() override;
   virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;
 
+  /**
+   * Compute the Jacobian of the Petrov-Galerkin stabilization addition with respect to the provided
+   * velocity component
+   */
+  Real computePGVelocityJacobian(unsigned short component);
+
+  /**
+   * Compute the strong residual, e.g. -rho * ddisp/dt * grad(u)
+   */
+  Real strongResidual();
+
   const VariableValue & _disp_x_dot;
   const VariableValue & _d_disp_x_dot;
   const unsigned int _disp_x_id;
@@ -37,4 +48,10 @@ class ConvectedMesh : public Kernel
   const unsigned int _disp_z_id;
 
   const MaterialProperty<Real> & _rho;
+
+  /// Whether this kernel should include Streamline-Upwind Petrov-Galerkin stabilization
+  const bool _supg;
+
+  /// The velocity component this kernel is acting on
+  unsigned short _component;
 };
diff --git a/modules/fsi/include/kernels/ConvectedMeshPSPG.h b/modules/fsi/include/kernels/ConvectedMeshPSPG.h
new file mode 100644
index 000000000000..617caeeb54cb
--- /dev/null
+++ b/modules/fsi/include/kernels/ConvectedMeshPSPG.h
@@ -0,0 +1,48 @@
+//* 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 "INSBase.h"
+
+/**
+ *  Computes residual and Jacobian contributions for the PSPG stabilization term for mesh advection
+ **/
+class ConvectedMeshPSPG : public INSBase
+{
+public:
+  static InputParameters validParams();
+
+  ConvectedMeshPSPG(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpResidual() override;
+  virtual Real computeQpJacobian() override;
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;
+
+  RealVectorValue dStrongResidualDDisp(unsigned short component);
+  RealVectorValue dStrongResidualDVel(unsigned short component);
+
+  /**
+   * Compute the strong residual, e.g. -rho * ddisp/dt * grad(u)
+   */
+  RealVectorValue strongResidual();
+
+  const VariableValue & _disp_x_dot;
+  const VariableValue & _d_disp_x_dot;
+  const unsigned int _disp_x_id;
+  const VariableValue & _disp_y_dot;
+  const VariableValue & _d_disp_y_dot;
+  const unsigned int _disp_y_id;
+  const VariableValue & _disp_z_dot;
+  const VariableValue & _d_disp_z_dot;
+  const unsigned int _disp_z_id;
+
+  const MaterialProperty<Real> & _rho;
+};
diff --git a/modules/fsi/src/interfacekernels/ADPenaltyVelocityContinuity.C b/modules/fsi/src/interfacekernels/ADPenaltyVelocityContinuity.C
new file mode 100644
index 000000000000..bd44f5906126
--- /dev/null
+++ b/modules/fsi/src/interfacekernels/ADPenaltyVelocityContinuity.C
@@ -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
+
+#include "ADPenaltyVelocityContinuity.h"
+#include "Assembly.h"
+
+registerMooseObject("FsiApp", ADPenaltyVelocityContinuity);
+
+InputParameters
+ADPenaltyVelocityContinuity::validParams()
+{
+  InputParameters params = InterfaceKernelBase::validParams();
+  params.addClassDescription(
+      "Enforces continuity of flux and continuity of solution via penalty across an interface.");
+  params.addRequiredParam<Real>(
+      "penalty",
+      "The penalty that penalizes jump between primary and neighbor secondary variables.");
+  params.addRequiredCoupledVar("fluid_velocity", "The fluid velocity");
+  params.renameCoupledVar("neighbor_var", "displacements", "All the displacement variables");
+  params.addRequiredCoupledVar("solid_velocities", "The solid velocity components");
+  return params;
+}
+
+ADPenaltyVelocityContinuity::ADPenaltyVelocityContinuity(const InputParameters & parameters)
+  : InterfaceKernelBase(parameters),
+    _penalty(getParam<Real>("penalty")),
+    _velocity_var(getVectorVar("fluid_velocity", 0)),
+    _velocity(adCoupledVectorValue("fluid_velocity")),
+    _ad_JxW(_assembly.adJxWFace()),
+    _ad_coord(_assembly.adCoordTransformation())
+{
+  if (isCoupledConstant("fluid_velocity"))
+    paramError("fluid_velocity", "The fluid velocity must be an actual variable");
+
+  _solid_velocities.resize(coupledComponents("solid_velocities"));
+  for (const auto i : index_range(_solid_velocities))
+    _solid_velocities[i] = &adCoupledNeighborValue("solid_velocities", i);
+
+  _displacements.resize(coupledComponents("displacements"));
+  for (const auto i : index_range(_displacements))
+    _displacements[i] = getVar("displacements", i);
+}
+
+void
+ADPenaltyVelocityContinuity::computeResidual()
+{
+  _qp_jumps.resize(_qrule->n_points());
+  for (auto & qp_jump : _qp_jumps)
+    qp_jump = 0;
+
+  const auto solid_velocity = [&](const auto qp)
+  {
+    ADRealVectorValue ret;
+    for (const auto i : index_range(_solid_velocities))
+      if (_solid_velocities[i])
+        ret(i) = (*_solid_velocities[i])[qp];
+
+    return ret;
+  };
+
+  for (const auto qp : make_range(_qrule->n_points()))
+    _qp_jumps[qp] = _ad_JxW[qp] * _ad_coord[qp] * _penalty * (_velocity[qp] - solid_velocity(qp));
+
+  // Fluid velocity residuals
+  {
+    const auto & phi = _velocity_var->phiFace();
+    const auto & dof_indices = _velocity_var->dofIndices();
+    mooseAssert(phi.size() == dof_indices.size(), "These should be the same");
+
+    _residuals.resize(phi.size());
+    for (auto & r : _residuals)
+      r = 0;
+
+    for (const auto i : index_range(phi))
+      for (const auto qp : make_range(_qrule->n_points()))
+        _residuals[i] += phi[i][qp] * _qp_jumps[qp];
+
+    addResidualsAndJacobian(_assembly, _residuals, dof_indices, _velocity_var->scalingFactor());
+  }
+
+  // Displacement residuals
+  for (const auto dim : index_range(_displacements))
+  {
+    const auto * const disp_var = _displacements[dim];
+    if (disp_var)
+    {
+      const auto & phi = disp_var->phiFaceNeighbor();
+      const auto & dof_indices = disp_var->dofIndicesNeighbor();
+      mooseAssert(phi.size() == dof_indices.size(), "These should be the same");
+
+      _residuals.resize(phi.size());
+      for (auto & r : _residuals)
+        r = 0;
+
+      for (const auto qp : make_range(_qrule->n_points()))
+        for (const auto i : index_range(phi))
+          _residuals[i] -= phi[i][qp] * _qp_jumps[qp](dim);
+
+      addResidualsAndJacobian(_assembly, _residuals, dof_indices, disp_var->scalingFactor());
+    }
+  }
+}
diff --git a/modules/fsi/src/kernels/ConvectedMesh.C b/modules/fsi/src/kernels/ConvectedMesh.C
index da66fe028942..43832bdc951c 100644
--- a/modules/fsi/src/kernels/ConvectedMesh.C
+++ b/modules/fsi/src/kernels/ConvectedMesh.C
@@ -14,18 +14,20 @@ registerMooseObject("FsiApp", ConvectedMesh);
 InputParameters
 ConvectedMesh::validParams()
 {
-  InputParameters params = Kernel::validParams();
+  InputParameters params = INSBase::validParams();
   params.addClassDescription(
       "Corrects the convective derivative for situations in which the fluid mesh is dynamic.");
   params.addRequiredCoupledVar("disp_x", "The x displacement");
   params.addCoupledVar("disp_y", "The y displacement");
   params.addCoupledVar("disp_z", "The z displacement");
   params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
+  params.addParam<bool>(
+      "supg", false, "Whether to perform SUPG stabilization of the momentum residuals");
   return params;
 }
 
 ConvectedMesh::ConvectedMesh(const InputParameters & parameters)
-  : Kernel(parameters),
+  : INSBase(parameters),
     _disp_x_dot(coupledDot("disp_x")),
     _d_disp_x_dot(coupledDotDu("disp_x")),
     _disp_x_id(coupled("disp_x")),
@@ -35,33 +37,97 @@ ConvectedMesh::ConvectedMesh(const InputParameters & parameters)
     _disp_z_dot(isCoupled("disp_z") ? coupledDot("disp_z") : _zero),
     _d_disp_z_dot(isCoupled("disp_z") ? coupledDotDu("disp_z") : _zero),
     _disp_z_id(coupled("disp_z")),
-    _rho(getMaterialProperty<Real>("rho_name"))
+    _rho(getMaterialProperty<Real>("rho_name")),
+    _supg(getParam<bool>("supg"))
 {
+  if (_var.number() == _u_vel_var_number)
+    _component = 0;
+  else if (_var.number() == _v_vel_var_number)
+    _component = 1;
+  else if (_var.number() == _w_vel_var_number)
+    _component = 2;
+  else
+    paramError("variable", "The variable must match one of the velocity variables.");
+}
+
+Real
+ConvectedMesh::strongResidual()
+{
+  return -_rho[_qp] * RealVectorValue(_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]) *
+         _grad_u[_qp];
 }
 
 Real
 ConvectedMesh::computeQpResidual()
 {
-  return _test[_i][_qp] * -_rho[_qp] *
-         RealVectorValue(_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]) * _grad_u[_qp];
+  auto test = _test[_i][_qp];
+  const auto U = relativeVelocity();
+  if (_supg)
+    test += tau() * _grad_test[_i][_qp] * U;
+  return test * strongResidual();
+}
+
+Real
+ConvectedMesh::computePGVelocityJacobian(const unsigned short component)
+{
+  const auto U = relativeVelocity();
+  return strongResidual() * ((dTauDUComp(component) * _grad_test[_i][_qp] * U) +
+                             (tau() * _grad_test[_i][_qp](component) * _phi[_j][_qp]));
 }
 
 Real
 ConvectedMesh::computeQpJacobian()
 {
-  return _test[_i][_qp] * -_rho[_qp] *
-         RealVectorValue(_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]) * _grad_phi[_j][_qp];
+  auto test = _test[_i][_qp];
+  const auto U = relativeVelocity();
+  if (_supg)
+    test += tau() * _grad_test[_i][_qp] * U;
+  auto jac = test * -_rho[_qp] *
+             RealVectorValue(_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]) *
+             _grad_phi[_j][_qp];
+  if (_supg)
+    jac += computePGVelocityJacobian(_component);
+
+  return jac;
 }
 
 Real
 ConvectedMesh::computeQpOffDiagJacobian(unsigned int jvar)
 {
+  mooseAssert(jvar != _var.number(), "Making sure I understand how old hand-coded Jacobians work.");
+
+  auto test = _test[_i][_qp];
+  const auto U = relativeVelocity();
+  if (_supg)
+    test += tau() * _grad_test[_i][_qp] * U;
+
   if (jvar == _disp_x_id)
-    return _test[_i][_qp] * -_rho[_qp] * _phi[_j][_qp] * _d_disp_x_dot[_qp] * _grad_u[_qp](0);
+    return test * -_rho[_qp] * _phi[_j][_qp] * _d_disp_x_dot[_qp] * _grad_u[_qp](0);
   else if (jvar == _disp_y_id)
-    return _test[_i][_qp] * -_rho[_qp] * _phi[_j][_qp] * _d_disp_y_dot[_qp] * _grad_u[_qp](1);
+    return test * -_rho[_qp] * _phi[_j][_qp] * _d_disp_y_dot[_qp] * _grad_u[_qp](1);
   else if (jvar == _disp_z_id)
-    return _test[_i][_qp] * -_rho[_qp] * _phi[_j][_qp] * _d_disp_z_dot[_qp] * _grad_u[_qp](2);
+    return test * -_rho[_qp] * _phi[_j][_qp] * _d_disp_z_dot[_qp] * _grad_u[_qp](2);
+  else if (jvar == _u_vel_var_number)
+  {
+    if (_supg)
+      return computePGVelocityJacobian(0);
+    else
+      return 0;
+  }
+  else if (jvar == _v_vel_var_number)
+  {
+    if (_supg)
+      return computePGVelocityJacobian(1);
+    else
+      return 0;
+  }
+  else if (jvar == _w_vel_var_number)
+  {
+    if (_supg)
+      return computePGVelocityJacobian(2);
+    else
+      return 0;
+  }
   else
     return 0.0;
 }
diff --git a/modules/fsi/src/kernels/ConvectedMeshPSPG.C b/modules/fsi/src/kernels/ConvectedMeshPSPG.C
new file mode 100644
index 000000000000..e8da5f120ae3
--- /dev/null
+++ b/modules/fsi/src/kernels/ConvectedMeshPSPG.C
@@ -0,0 +1,127 @@
+//* 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 "ConvectedMeshPSPG.h"
+
+registerMooseObject("FsiApp", ConvectedMeshPSPG);
+
+InputParameters
+ConvectedMeshPSPG::validParams()
+{
+  InputParameters params = INSBase::validParams();
+  params.addClassDescription(
+      "Corrects the convective derivative for situations in which the fluid mesh is dynamic.");
+  params.addRequiredCoupledVar("disp_x", "The x displacement");
+  params.addCoupledVar("disp_y", "The y displacement");
+  params.addCoupledVar("disp_z", "The z displacement");
+  params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
+  return params;
+}
+
+ConvectedMeshPSPG::ConvectedMeshPSPG(const InputParameters & parameters)
+  : INSBase(parameters),
+    _disp_x_dot(coupledDot("disp_x")),
+    _d_disp_x_dot(coupledDotDu("disp_x")),
+    _disp_x_id(coupled("disp_x")),
+    _disp_y_dot(isCoupled("disp_y") ? coupledDot("disp_y") : _zero),
+    _d_disp_y_dot(isCoupled("disp_y") ? coupledDotDu("disp_y") : _zero),
+    _disp_y_id(coupled("disp_y")),
+    _disp_z_dot(isCoupled("disp_z") ? coupledDot("disp_z") : _zero),
+    _d_disp_z_dot(isCoupled("disp_z") ? coupledDotDu("disp_z") : _zero),
+    _disp_z_id(coupled("disp_z")),
+    _rho(getMaterialProperty<Real>("rho_name"))
+{
+}
+
+RealVectorValue
+ConvectedMeshPSPG::strongResidual()
+{
+  const auto minus_rho_ddisp_dt =
+      -_rho[_qp] * RealVectorValue(_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]);
+  return RealVectorValue(minus_rho_ddisp_dt * _grad_u_vel[_qp],
+                         minus_rho_ddisp_dt * _grad_v_vel[_qp],
+                         minus_rho_ddisp_dt * _grad_w_vel[_qp]);
+}
+
+RealVectorValue
+ConvectedMeshPSPG::dStrongResidualDDisp(const unsigned short component)
+{
+  const auto & ddisp_dot = [&]() -> const VariableValue &
+  {
+    switch (component)
+    {
+      case 0:
+        return _d_disp_x_dot;
+      case 1:
+        return _d_disp_y_dot;
+      case 2:
+        return _d_disp_z_dot;
+      default:
+        mooseError("Invalid component");
+    }
+  }();
+
+  // Only non-zero component will be from 'component'
+  RealVectorValue ddisp_dt;
+  ddisp_dt(component) = _phi[_j][_qp] * ddisp_dot[_qp];
+
+  const auto minus_rho_ddisp_dt = -_rho[_qp] * ddisp_dt;
+  return RealVectorValue(minus_rho_ddisp_dt * _grad_u_vel[_qp],
+                         minus_rho_ddisp_dt * _grad_v_vel[_qp],
+                         minus_rho_ddisp_dt * _grad_w_vel[_qp]);
+}
+
+RealVectorValue
+ConvectedMeshPSPG::dStrongResidualDVel(const unsigned short component)
+{
+  const auto minus_rho_ddisp_dt =
+      -_rho[_qp] * RealVectorValue(_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]);
+
+  // Only non-zero component will be from 'component'
+  RealVectorValue ret;
+  ret(component) = minus_rho_ddisp_dt * _grad_phi[_j][_qp];
+  return ret;
+}
+
+Real
+ConvectedMeshPSPG::computeQpResidual()
+{
+  return -tau() / _rho[_qp] * _grad_test[_i][_qp] * strongResidual();
+}
+
+Real
+ConvectedMeshPSPG::computeQpJacobian()
+{
+  // No derivative with respect to pressure
+  return 0;
+}
+
+Real
+ConvectedMeshPSPG::computeQpOffDiagJacobian(unsigned int jvar)
+{
+  mooseAssert(jvar != _var.number(), "Making sure I understand how old hand-coded Jacobians work.");
+
+  if (jvar == _disp_x_id)
+    return -tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDDisp(0);
+  else if (jvar == _disp_y_id)
+    return -tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDDisp(1);
+  else if (jvar == _disp_z_id)
+    return -tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDDisp(2);
+  else if (jvar == _u_vel_var_number)
+    return -dTauDUComp(0) / _rho[_qp] * _grad_test[_i][_qp] * strongResidual() -
+           tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDVel(0);
+  else if (jvar == _v_vel_var_number)
+    return -dTauDUComp(1) / _rho[_qp] * _grad_test[_i][_qp] * strongResidual() -
+           tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDVel(1);
+  else if (jvar == _w_vel_var_number)
+    return -dTauDUComp(2) / _rho[_qp] * _grad_test[_i][_qp] * strongResidual() -
+           tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDVel(2);
+  else
+    return 0.0;
+}
diff --git a/modules/fsi/test/tests/2d-small-strain-transient/ad-fsi-flat-channel.i b/modules/fsi/test/tests/2d-small-strain-transient/ad-fsi-flat-channel.i
new file mode 100644
index 000000000000..a605cbb299fe
--- /dev/null
+++ b/modules/fsi/test/tests/2d-small-strain-transient/ad-fsi-flat-channel.i
@@ -0,0 +1,284 @@
+[GlobalParams]
+  displacements = 'disp_x disp_y'
+  order = FIRST
+  preset = false
+  use_displaced_mesh = true
+[]
+
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 2
+    xmin = 0
+    xmax = 3.0
+    ymin = 0
+    ymax = 1.0
+    nx = 10
+    ny = 15
+    elem_type = QUAD4
+  []
+  [subdomain1]
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '0.0 0.5 0'
+    block_id = 1
+    top_right = '3.0 1.0 0'
+    input = gmg
+  []
+  [interface]
+    type = SideSetsBetweenSubdomainsGenerator
+    primary_block = '0'
+    paired_block = '1'
+    new_boundary = 'master0_interface'
+    input = subdomain1
+  []
+  [break_boundary]
+    type = BreakBoundaryOnSubdomainGenerator
+    input = interface
+  []
+[]
+
+[Variables]
+  [vel]
+    block = 0
+    family = LAGRANGE_VEC
+  []
+  [p]
+    block = 0
+    order = FIRST
+  []
+  [disp_x]
+  []
+  [disp_y]
+  []
+  [vel_x_solid]
+    block = 1
+  []
+  [vel_y_solid]
+    block = 1
+  []
+[]
+
+[Kernels]
+  [mass]
+    type = INSADMass
+    variable = p
+    block = 0
+  []
+  [mass_pspg]
+    type = INSADMassPSPG
+    variable = p
+    block = 0
+  []
+  [momentum_time]
+    type = INSADMomentumTimeDerivative
+    variable = vel
+    block = 0
+  []
+  [momentum_convection]
+    type = INSADMomentumAdvection
+    variable = vel
+    block = 0
+  []
+  [momentum_viscous]
+    type = INSADMomentumViscous
+    variable = vel
+    block = 0
+  []
+  [momentum_pressure]
+    type = INSADMomentumPressure
+    variable = vel
+    pressure = p
+    integrate_p_by_parts = true
+    block = 0
+  []
+  [momentum_supg]
+    type = INSADMomentumSUPG
+    variable = vel
+    material_velocity = relative_velocity
+    block = 0
+  []
+  [momentum_mesh]
+    type = INSADMomentumMeshAdvection
+    variable = vel
+    disp_x = 'disp_x'
+    disp_y = 'disp_y'
+    block = 0
+  []
+  [disp_x_fluid]
+    type = Diffusion
+    variable = disp_x
+    block = 0
+    use_displaced_mesh = false
+  []
+  [disp_y_fluid]
+    type = Diffusion
+    variable = disp_y
+    block = 0
+    use_displaced_mesh = false
+  []
+  [accel_tensor_x]
+    type = CoupledTimeDerivative
+    variable = disp_x
+    v = vel_x_solid
+    block = 1
+    use_displaced_mesh = false
+  []
+  [accel_tensor_y]
+    type = CoupledTimeDerivative
+    variable = disp_y
+    v = vel_y_solid
+    block = 1
+    use_displaced_mesh = false
+  []
+  [vxs_time_derivative_term]
+    type = CoupledTimeDerivative
+    variable = vel_x_solid
+    v = disp_x
+    block = 1
+    use_displaced_mesh = false
+  []
+  [vys_time_derivative_term]
+    type = CoupledTimeDerivative
+    variable = vel_y_solid
+    v = disp_y
+    block = 1
+    use_displaced_mesh = false
+  []
+  [source_vxs]
+    type = MatReaction
+    variable = vel_x_solid
+    block = 1
+    mob_name = 1
+    use_displaced_mesh = false
+  []
+  [source_vys]
+    type = MatReaction
+    variable = vel_y_solid
+    block = 1
+    mob_name = 1
+    use_displaced_mesh = false
+  []
+[]
+
+[InterfaceKernels]
+  [penalty]
+    type = ADPenaltyVelocityContinuity
+    variable = vel
+    fluid_velocity = vel
+    displacements = 'disp_x disp_y'
+    solid_velocities = 'vel_x_solid vel_y_solid'
+    boundary = master0_interface
+    penalty = 1e6
+  []
+[]
+
+[Modules/TensorMechanics/Master]
+  [solid_domain]
+    strain = SMALL
+    incremental = false
+    # generate_output = 'strain_xx strain_yy strain_zz' ## Not at all necessary, but nice
+    block = '1'
+  []
+[]
+
+[Materials]
+  [elasticity_tensor]
+    type = ComputeIsotropicElasticityTensor
+    youngs_modulus = 1e2
+    poissons_ratio = 0.3
+    block = '1'
+    use_displaced_mesh = false
+  []
+  [small_stress]
+    type = ComputeLinearElasticStress
+    block = 1
+  []
+  [const]
+    type = ADGenericConstantMaterial
+    block = 0
+    prop_names = 'rho mu'
+    prop_values = '1  1'
+  []
+  [ins_mat]
+    type = INSADTauMaterial
+    velocity = vel
+    pressure = p
+    block = 0
+  []
+[]
+
+[BCs]
+  [fluid_bottom]
+    type = ADVectorFunctionDirichletBC
+    variable = vel
+    boundary = 'bottom'
+    function_x = 0
+    function_y = 0
+  []
+  [fluid_left]
+    type = ADVectorFunctionDirichletBC
+    variable = vel
+    boundary = 'left_to_0'
+    function_x = 'inlet_func'
+    function_y = 0
+    # The displacements actually affect the result of the function evaluation so in order to eliminate the impact
+    # on the Jacobian we set 'use_displaced_mesh = false' here
+    use_displaced_mesh = false
+  []
+  [no_disp_x]
+    type = DirichletBC
+    variable = disp_x
+    boundary = 'bottom top left_to_1 right_to_1 left_to_0 right_to_0'
+    value = 0
+  []
+  [no_disp_y]
+    type = DirichletBC
+    variable = disp_y
+    boundary = 'bottom top left_to_1 right_to_1 left_to_0 right_to_0'
+    value = 0
+  []
+  [solid_x_no_slip]
+    type = DirichletBC
+    variable = vel_x_solid
+    boundary = 'top left_to_1 right_to_1'
+    value = 0.0
+  []
+  [solid_y_no_slip]
+    type = DirichletBC
+    variable = vel_y_solid
+    boundary = 'top left_to_1 right_to_1'
+    value = 0.0
+  []
+[]
+
+[Preconditioning]
+  [SMP]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = 5
+  # num_steps = 60
+  dt = 0.1
+  dtmin = 0.1
+  solve_type = 'NEWTON'
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu       NONZERO'
+  line_search = none
+  nl_rel_tol = 1e-50
+  nl_abs_tol = 1e-10
+[]
+
+[Outputs]
+  exodus = true
+[]
+
+[Functions]
+  [inlet_func]
+    type = ParsedFunction
+    expression = '(-16 * (y - 0.25)^2 + 1) * (1 + cos(t))'
+  []
+[]
diff --git a/modules/fsi/test/tests/2d-small-strain-transient/custom.cmp b/modules/fsi/test/tests/2d-small-strain-transient/custom.cmp
new file mode 100644
index 000000000000..5933c52784c2
--- /dev/null
+++ b/modules/fsi/test/tests/2d-small-strain-transient/custom.cmp
@@ -0,0 +1,43 @@
+
+#  *****************************************************************
+#             EXODIFF	(Version: 2.90) Modified: 2018-02-15
+#             Authors:  Richard Drake, rrdrake@sandia.gov
+#                       Greg Sjaardema, gdsjaar@sandia.gov
+#             Run on    2023/04/26   11:42:41 PDT
+#  *****************************************************************
+
+#  FILE 1: /home/lindad/projects/moose/modules/fsi/test/tests/2d-small-strain-transient/ad-fsi-flat-channel_out.e
+#   Title: ad-fsi-flat-channel_out.e
+#          Dim = 2, Blocks = 2, Nodes = 35, Elements = 6, Nodesets = 12, Sidesets = 12
+#          Vars: Global = 0, Nodal = 7, Element = 0, Nodeset = 0, Sideset = 0, Times = 2
+
+
+# ==============================================================
+#  NOTE: All node and element ids are reported as global ids.
+
+# NOTES:  - The min/max values are reporting the min/max in absolute value.
+#         - Time values (t) are 1-offset time step numbers.
+#         - Element block numbers are the block ids.
+#         - Node(n) and element(e) numbers are 1-offset.
+
+COORDINATES absolute 1.e-6    # min separation not calculated
+
+TIME STEPS relative 1.e-6 floor 0.0     # min:               0 @ t1 max:             0.1 @ t2
+
+
+# No GLOBAL VARIABLES
+
+NODAL VARIABLES relative 5.5e-6 floor 1e-10
+	disp_x       # min:               0 @ t1,n1	max:     0.035855951 @ t2,n7
+	disp_y       # min:               0 @ t1,n1	max:     0.065177301 @ t2,n7
+	p            # min:               0 @ t1,n1	max:       38.023022 @ t2,n4
+	vel_x        # min:               0 @ t1,n1	max:        1.773337 @ t2,n4
+	vel_x_solid  # min:               0 @ t1,n1	max:      0.35855951 @ t2,n7
+	vel_y        # min:               0 @ t1,n1	max:         0.65179 @ t2,n7
+	vel_y_solid  # min:               0 @ t1,n1	max:      0.65177301 @ t2,n7
+
+# No ELEMENT VARIABLES
+
+# No NODESET VARIABLES
+
+# No SIDESET VARIABLES
diff --git a/modules/fsi/test/tests/2d-small-strain-transient/fsi_flat_channel.i b/modules/fsi/test/tests/2d-small-strain-transient/fsi_flat_channel.i
index 41e9b6875613..81b65334510f 100644
--- a/modules/fsi/test/tests/2d-small-strain-transient/fsi_flat_channel.i
+++ b/modules/fsi/test/tests/2d-small-strain-transient/fsi_flat_channel.i
@@ -5,20 +5,24 @@
   convective_term = true
   transient_term = true
   pspg = true
+  supg = true
   displacements = 'disp_x disp_y'
+  preset = false
+  order = FIRST
+  use_displaced_mesh = true
 []
 
 [Mesh]
   [gmg]
-  type = GeneratedMeshGenerator
-  dim = 2
-  xmin = 0
-  xmax = 3.0
-  ymin = 0
-  ymax = 1.0
-  nx = 10
-  ny = 15
-  elem_type = QUAD4
+    type = GeneratedMeshGenerator
+    dim = 2
+    xmin = 0
+    xmax = 3.0
+    ymin = 0
+    ymax = 1.0
+    nx = 10
+    ny = 15
+    elem_type = QUAD4
   []
 
   [subdomain1]
@@ -70,13 +74,11 @@
     type = INSMomentumTimeDerivative
     variable = vel_x
     block = 0
-    use_displaced_mesh = true
   [../]
   [./vel_y_time]
     type = INSMomentumTimeDerivative
     variable = vel_y
     block = 0
-    use_displaced_mesh = true
   [../]
   [./mass]
     type = INSMass
@@ -85,7 +87,8 @@
     v = vel_y
     pressure = p
     block = 0
-    use_displaced_mesh = true
+    disp_x = disp_x
+    disp_y = disp_y
   [../]
   [./x_momentum_space]
     type = INSMomentumLaplaceForm
@@ -95,7 +98,8 @@
     pressure = p
     component = 0
     block = 0
-    use_displaced_mesh = true
+    disp_x = disp_x
+    disp_y = disp_y
   [../]
   [./y_momentum_space]
     type = INSMomentumLaplaceForm
@@ -105,69 +109,92 @@
     pressure = p
     component = 1
     block = 0
-    use_displaced_mesh = true
+    disp_x = disp_x
+    disp_y = disp_y
   [../]
   [./vel_x_mesh]
     type = ConvectedMesh
     disp_x = disp_x
     disp_y = disp_y
     variable = vel_x
+    u = vel_x
+    v = vel_y
+    pressure = p
     block = 0
-    use_displaced_mesh = true
   [../]
   [./vel_y_mesh]
     type = ConvectedMesh
     disp_x = disp_x
     disp_y = disp_y
     variable = vel_y
+    u = vel_x
+    v = vel_y
+    pressure = p
+    block = 0
+  [../]
+  [./p_mesh]
+    type = ConvectedMeshPSPG
+    disp_x = disp_x
+    disp_y = disp_y
+    variable = p
+    u = vel_x
+    v = vel_y
+    pressure = p
     block = 0
-    use_displaced_mesh = true
   [../]
   [./disp_x_fluid]
     type = Diffusion
     variable = disp_x
     block = 0
+    use_displaced_mesh = false
   [../]
   [./disp_y_fluid]
     type = Diffusion
     variable = disp_y
     block = 0
+    use_displaced_mesh = false
   [../]
   [./accel_tensor_x]
     type = CoupledTimeDerivative
     variable = disp_x
     v = vel_x_solid
     block = 1
+    use_displaced_mesh = false
   [../]
   [./accel_tensor_y]
     type = CoupledTimeDerivative
     variable = disp_y
     v = vel_y_solid
     block = 1
+    use_displaced_mesh = false
   [../]
   [./vxs_time_derivative_term]
     type = CoupledTimeDerivative
     variable = vel_x_solid
     v = disp_x
     block = 1
+    use_displaced_mesh = false
   [../]
   [./vys_time_derivative_term]
     type = CoupledTimeDerivative
     variable = vel_y_solid
     v = disp_y
     block = 1
+    use_displaced_mesh = false
   [../]
   [./source_vxs]
     type = MatReaction
     variable = vel_x_solid
     block = 1
     mob_name = 1
+    use_displaced_mesh = false
   [../]
   [./source_vys]
     type = MatReaction
     variable = vel_y_solid
     block = 1
     mob_name = 1
+    use_displaced_mesh = false
   [../]
 []
 
@@ -205,6 +232,7 @@
     youngs_modulus = 1e2
     poissons_ratio = 0.3
     block = '1'
+    use_displaced_mesh = false
   [../]
   [./small_stress]
     type = ComputeLinearElasticStress
@@ -215,6 +243,7 @@
     block = 0
     prop_names = 'rho mu'
     prop_values = '1  1'
+    use_displaced_mesh = false
   [../]
 []
 
@@ -277,9 +306,11 @@
   dt = 0.1
   dtmin = 0.1
   solve_type = 'PJFNK'
-  petsc_options_iname = '-pc_type'
-  petsc_options_value = 'lu'
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu       NONZERO'
   line_search = none
+  nl_rel_tol = 1e-50
+  nl_abs_tol = 1e-10
 []
 
 [Outputs]
diff --git a/modules/fsi/test/tests/2d-small-strain-transient/gold/ad-fsi-flat-channel_out.e b/modules/fsi/test/tests/2d-small-strain-transient/gold/ad-fsi-flat-channel_out.e
new file mode 120000
index 000000000000..ffa472d109a3
--- /dev/null
+++ b/modules/fsi/test/tests/2d-small-strain-transient/gold/ad-fsi-flat-channel_out.e
@@ -0,0 +1 @@
+fsi_flat_channel_out.e
\ No newline at end of file
diff --git a/modules/fsi/test/tests/2d-small-strain-transient/gold/fsi_flat_channel_out.e b/modules/fsi/test/tests/2d-small-strain-transient/gold/fsi_flat_channel_out.e
index 0d72cc5b1c88..9678cb1f2d63 100644
Binary files a/modules/fsi/test/tests/2d-small-strain-transient/gold/fsi_flat_channel_out.e and b/modules/fsi/test/tests/2d-small-strain-transient/gold/fsi_flat_channel_out.e differ
diff --git a/modules/fsi/test/tests/2d-small-strain-transient/tests b/modules/fsi/test/tests/2d-small-strain-transient/tests
index 92074d61b10f..71037c3d3e4f 100644
--- a/modules/fsi/test/tests/2d-small-strain-transient/tests
+++ b/modules/fsi/test/tests/2d-small-strain-transient/tests
@@ -1,14 +1,42 @@
 [Tests]
   design = 'modules/fsi/index.md'
-
-  [./fsi]
-    type = 'Exodiff'
-    input = 'fsi_flat_channel.i'
-    exodiff = 'fsi_flat_channel_out.e'
-    requirement = 'We shall be able to combine navier-stokes and tensor-mechanics simulation capabilities to do some basic fluid-structure interaction simulations'
-    issues = '#12462'
-    # vel_x_solid is not accurate
-    # We mmay improve that in the future
-    rel_err = 1e-4
-  [../]
+  issues = '#12462'
+  [fsi]
+    requirement = 'We shall be able to combine navier-stokes and tensor-mechanics simulation capabilities to do some basic fluid-structure interaction simulations and produce equivalent results when the Navier-Stokes equations are implemented using'
+    [INS]
+      type = 'Exodiff'
+      input = 'fsi_flat_channel.i'
+      exodiff = 'fsi_flat_channel_out.e'
+      detail = 'scalar field variables for velocity components and hand-coded Jacobians,'
+      # vel_x_solid is not accurate
+      # We may improve that in the future
+      rel_err = 1e-4
+    []
+    [INSAD]
+      type = 'Exodiff'
+      input = 'ad-fsi-flat-channel.i'
+      exodiff = 'ad-fsi-flat-channel_out.e'
+      detail = 'and a vector field variable for velocity and a Jacobian computed using automatic differentation. Additionally,'
+      # vel_x_solid is not accurate
+      # We mmay improve that in the future
+      rel_err = 1e-4
+    []
+    [INSAD_displaced_jac]
+      type = 'PetscJacobianTester'
+      input = 'ad-fsi-flat-channel.i'
+      ratio_tol = 1e-6
+      difference_tol = 1
+      run_sim = True
+      cli_args = 'Mesh/gmg/nx=2 Mesh/gmg/ny=3 Executioner/num_steps=1'
+      detail = 'the automatic differentation Jacobian shall be (nearly) perfect when the fluid domain equations are run on the displaced mesh, and'
+    []
+    [INSAD_undisplaced_jac]
+      type = 'PetscJacobianTester'
+      input = 'ad-fsi-flat-channel.i'
+      difference_tol = 1
+      run_sim = True
+      cli_args = 'Mesh/gmg/nx=2 Mesh/gmg/ny=3 Executioner/num_steps=1 GlobalParams/use_displaced_mesh=false'
+      detail = 'perfect when the fluid domain equations are run on the undisplaced mesh, and'
+    []
+  []
 []
diff --git a/modules/navier_stokes/examples/laser-welding/3d.i b/modules/navier_stokes/examples/laser-welding/3d.i
index 2000fbb842d1..70a140f2496a 100644
--- a/modules/navier_stokes/examples/laser-welding/3d.i
+++ b/modules/navier_stokes/examples/laser-welding/3d.i
@@ -104,7 +104,7 @@ sb=5.67e-8
   [momentum_supg]
     type = INSADMomentumSUPG
     variable = vel
-    velocity = vel
+    material_velocity = relative_velocity
     use_displaced_mesh = true
   []
   [temperature_time]
diff --git a/modules/navier_stokes/include/kernels/INSADEnergyMeshAdvection.h b/modules/navier_stokes/include/kernels/INSADEnergyMeshAdvection.h
index 736fa1a9ed37..533fa314fbc1 100644
--- a/modules/navier_stokes/include/kernels/INSADEnergyMeshAdvection.h
+++ b/modules/navier_stokes/include/kernels/INSADEnergyMeshAdvection.h
@@ -25,7 +25,7 @@ class INSADEnergyMeshAdvection : public ADKernelValue
 protected:
   virtual ADReal precomputeQpResidual() override;
 
-  const ADMaterialProperty<Real> & _temperature_convected_mesh_strong_residual;
+  const ADMaterialProperty<Real> & _temperature_advected_mesh_strong_residual;
 
   friend class INSADMomentumMeshAdvection;
 };
diff --git a/modules/navier_stokes/include/kernels/INSADMomentumMeshAdvection.h b/modules/navier_stokes/include/kernels/INSADMomentumMeshAdvection.h
index 39e7d376260b..c42bac014416 100644
--- a/modules/navier_stokes/include/kernels/INSADMomentumMeshAdvection.h
+++ b/modules/navier_stokes/include/kernels/INSADMomentumMeshAdvection.h
@@ -30,7 +30,7 @@ class INSADMomentumMeshAdvection : public ADVectorKernelValue
 
   /// The strong residual for this object, computed by material classes to eliminate computation
   /// duplication in simulations in which we have stabilization
-  const ADMaterialProperty<RealVectorValue> & _convected_mesh_strong_residual;
+  const ADMaterialProperty<RealVectorValue> & _advected_mesh_strong_residual;
 };
 
 template <typename T>
@@ -57,7 +57,7 @@ INSADMomentumMeshAdvection::setDisplacementParams(T & mesh_convection_obj)
         "ins_ad_object_tracker");
     for (const auto block_id : mesh_convection_obj.blockIDs())
     {
-      obj_tracker.set("has_convected_mesh", true, block_id);
+      obj_tracker.set("has_advected_mesh", true, block_id);
       obj_tracker.set("disp_x", mesh_convection_obj.coupledName("disp_x"), block_id);
       if (mesh_convection_obj.isParamValid("disp_y"))
         obj_tracker.set("disp_y", mesh_convection_obj.coupledName("disp_y"), block_id);
diff --git a/modules/navier_stokes/include/kernels/INSBase.h b/modules/navier_stokes/include/kernels/INSBase.h
index 2efa3a85598a..1e34c800f52b 100644
--- a/modules/navier_stokes/include/kernels/INSBase.h
+++ b/modules/navier_stokes/include/kernels/INSBase.h
@@ -11,8 +11,6 @@
 
 #include "Kernel.h"
 
-// Forward Declarations
-
 /**
  * This class computes strong and weak components of the INS governing
  * equations.  These terms can then be assembled in child classes
@@ -60,6 +58,12 @@ class INSBase : public Kernel
   /// Provides tau which yields superconvergence for 1D advection-diffusion
   virtual Real tauNodal();
 
+  /**
+   * Compute the velocity. If displacements are provided, then this routine will subtract the time
+   * derivative of the displacements, e.g. the mesh velocity
+   */
+  RealVectorValue relativeVelocity() const;
+
   /// second derivatives of the shape function
   const VariablePhiSecond & _second_phi;
 
@@ -106,4 +110,13 @@ class INSBase : public Kernel
   bool _laplace;
   bool _convective_term;
   bool _transient_term;
+
+  /// Whether displacements are provided
+  const bool _disps_provided;
+  /// Time derivative of the x-displacement, mesh velocity in the x-direction
+  const VariableValue & _disp_x_dot;
+  /// Time derivative of the y-displacement, mesh velocity in the y-direction
+  const VariableValue & _disp_y_dot;
+  /// Time derivative of the z-displacement, mesh velocity in the z-direction
+  const VariableValue & _disp_z_dot;
 };
diff --git a/modules/navier_stokes/include/materials/INSAD3Eqn.h b/modules/navier_stokes/include/materials/INSAD3Eqn.h
index e29adea79233..2bd76f6461c3 100644
--- a/modules/navier_stokes/include/materials/INSAD3Eqn.h
+++ b/modules/navier_stokes/include/materials/INSAD3Eqn.h
@@ -39,7 +39,7 @@ class INSAD3Eqn : public INSADMaterial
 
   /// The strong residual for the temperature transport term corresponding to mesh velocity in an
   /// ALE simulation
-  ADMaterialProperty<Real> & _temperature_convected_mesh_strong_residual;
+  ADMaterialProperty<Real> & _temperature_advected_mesh_strong_residual;
 
   bool _has_ambient_convection;
   Real _ambient_convection_alpha;
diff --git a/modules/navier_stokes/include/materials/INSADMaterial.h b/modules/navier_stokes/include/materials/INSADMaterial.h
index f8f68f2255df..bc24d98add76 100644
--- a/modules/navier_stokes/include/materials/INSADMaterial.h
+++ b/modules/navier_stokes/include/materials/INSADMaterial.h
@@ -65,8 +65,11 @@ class INSADMaterial : public Material
   /// The mesh velocity
   ADMaterialProperty<RealVectorValue> & _mesh_velocity;
 
-  /// Strong residual corresponding to convected mesh term
-  ADMaterialProperty<RealVectorValue> & _convected_mesh_strong_residual;
+  /// The relative velocity, e.g. velocity - mesh_velocity
+  ADMaterialProperty<RealVectorValue> & _relative_velocity;
+
+  /// Strong residual corresponding to advected mesh term
+  ADMaterialProperty<RealVectorValue> & _advected_mesh_strong_residual;
 
   // /// Future addition pending addition of INSADMMSKernel.
   // /// Strong residual corresponding to the mms function term
@@ -124,7 +127,7 @@ class INSADMaterial : public Material
   std::vector<const Function *> _coupled_force_vector_function;
 
   /// Whether we have mesh convection
-  bool _has_convected_mesh;
+  bool _has_advected_mesh;
 
   /// The time derivative with respect to x-displacement
   const ADVariableValue * _disp_x_dot;
diff --git a/modules/navier_stokes/include/materials/INSADStabilized3Eqn.h b/modules/navier_stokes/include/materials/INSADStabilized3Eqn.h
index 338fd8de8bbe..359fdfeb1c69 100644
--- a/modules/navier_stokes/include/materials/INSADStabilized3Eqn.h
+++ b/modules/navier_stokes/include/materials/INSADStabilized3Eqn.h
@@ -9,19 +9,6 @@
 
 #pragma once
 
-#include "InputParameters.h"
-#include "NonlinearSystemBase.h"
-#include "FEProblemBase.h"
-#include "MaterialProperty.h"
-#include "MooseArray.h"
-
-#include "libmesh/elem.h"
-
-#include <vector>
-
-class INSADMaterial;
-class INSAD3Eqn;
-
 #include "INSADTauMaterial.h"
 #include "INSAD3Eqn.h"
 
@@ -42,13 +29,4 @@ class INSADStabilized3Eqn : public INSADTauMaterialTempl<INSAD3Eqn>
 
   ADMaterialProperty<Real> & _tau_energy;
   ADMaterialProperty<Real> & _temperature_strong_residual;
-
-  using INSADTauMaterialTempl<INSAD3Eqn>::_cp;
-  using INSADTauMaterialTempl<INSAD3Eqn>::_temperature_advective_strong_residual;
-  using INSADTauMaterialTempl<INSAD3Eqn>::_temperature_td_strong_residual;
-  using INSADTauMaterialTempl<INSAD3Eqn>::_temperature_ambient_convection_strong_residual;
-  using INSADTauMaterialTempl<INSAD3Eqn>::_temperature_source_strong_residual;
-  using INSADTauMaterialTempl<INSAD3Eqn>::_has_ambient_convection;
-  using INSADTauMaterialTempl<INSAD3Eqn>::_has_heat_source;
-  using INSADTauMaterialTempl<INSAD3Eqn>::_has_energy_transient;
 };
diff --git a/modules/navier_stokes/include/materials/INSADTauMaterial.h b/modules/navier_stokes/include/materials/INSADTauMaterial.h
index 9f6f1e8f2255..78994614ee87 100644
--- a/modules/navier_stokes/include/materials/INSADTauMaterial.h
+++ b/modules/navier_stokes/include/materials/INSADTauMaterial.h
@@ -90,11 +90,15 @@ class INSADTauMaterialTempl : public T
   /// The velocity system number
   const unsigned int _vel_sys_number;
 
+  /// The speed of the medium. This is the norm of the relative velocity, e.g. the velocity minus
+  /// the mesh velocity, at the current _qp
+  ADReal _speed;
+
   using T::_ad_q_point;
+  using T::_advected_mesh_strong_residual;
   using T::_advective_strong_residual;
   using T::_assembly;
   using T::_boussinesq_strong_residual;
-  using T::_convected_mesh_strong_residual;
   using T::_coord_sys;
   using T::_coupled_force_strong_residual;
   using T::_current_elem;
@@ -109,8 +113,8 @@ class INSADTauMaterialTempl : public T
   using T::_grad_p;
   using T::_grad_velocity;
   using T::_gravity_strong_residual;
+  using T::_has_advected_mesh;
   using T::_has_boussinesq;
-  using T::_has_convected_mesh;
   using T::_has_coupled_force;
   using T::_has_gravity;
   using T::_has_transient;
@@ -120,6 +124,7 @@ class INSADTauMaterialTempl : public T
   using T::_q_point;
   using T::_qp;
   using T::_qrule;
+  using T::_relative_velocity;
   using T::_rho;
   using T::_rz_axial_coord;
   using T::_rz_radial_coord;
@@ -335,8 +340,8 @@ INSADTauMaterialTempl<T>::computeQpProperties()
 
   const auto nu = _mu[_qp] / _rho[_qp];
   const auto transient_part = _has_transient ? 4. / (_dt * _dt) : 0.;
-  const auto speed = NS::computeSpeed(_velocity[_qp]);
-  _tau[_qp] = _alpha / std::sqrt(transient_part + (2. * speed / _hmax) * (2. * speed / _hmax) +
+  _speed = NS::computeSpeed(_relative_velocity[_qp]);
+  _tau[_qp] = _alpha / std::sqrt(transient_part + (2. * _speed / _hmax) * (2. * _speed / _hmax) +
                                  9. * (4. * nu / (_hmax * _hmax)) * (4. * nu / (_hmax * _hmax)));
 
   _momentum_strong_residual[_qp] =
@@ -351,8 +356,8 @@ INSADTauMaterialTempl<T>::computeQpProperties()
   if (_has_boussinesq)
     _momentum_strong_residual[_qp] += _boussinesq_strong_residual[_qp];
 
-  if (_has_convected_mesh)
-    _momentum_strong_residual[_qp] += _convected_mesh_strong_residual[_qp];
+  if (_has_advected_mesh)
+    _momentum_strong_residual[_qp] += _advected_mesh_strong_residual[_qp];
 
   if (_has_coupled_force)
     _momentum_strong_residual[_qp] += _coupled_force_strong_residual[_qp];
diff --git a/modules/navier_stokes/src/kernels/INSADEnergyMeshAdvection.C b/modules/navier_stokes/src/kernels/INSADEnergyMeshAdvection.C
index 3e87e1c1c56d..66764089f25f 100644
--- a/modules/navier_stokes/src/kernels/INSADEnergyMeshAdvection.C
+++ b/modules/navier_stokes/src/kernels/INSADEnergyMeshAdvection.C
@@ -24,8 +24,8 @@ INSADEnergyMeshAdvection::validParams()
 
 INSADEnergyMeshAdvection::INSADEnergyMeshAdvection(const InputParameters & parameters)
   : ADKernelValue(parameters),
-    _temperature_convected_mesh_strong_residual(
-        getADMaterialProperty<Real>("temperature_convected_mesh_strong_residual"))
+    _temperature_advected_mesh_strong_residual(
+        getADMaterialProperty<Real>("temperature_advected_mesh_strong_residual"))
 {
   INSADMomentumMeshAdvection::setDisplacementParams(*this);
 }
@@ -33,5 +33,5 @@ INSADEnergyMeshAdvection::INSADEnergyMeshAdvection(const InputParameters & param
 ADReal
 INSADEnergyMeshAdvection::precomputeQpResidual()
 {
-  return _temperature_convected_mesh_strong_residual[_qp];
+  return _temperature_advected_mesh_strong_residual[_qp];
 }
diff --git a/modules/navier_stokes/src/kernels/INSADMomentumMeshAdvection.C b/modules/navier_stokes/src/kernels/INSADMomentumMeshAdvection.C
index 2173720f247b..bcdbe0d12f4d 100644
--- a/modules/navier_stokes/src/kernels/INSADMomentumMeshAdvection.C
+++ b/modules/navier_stokes/src/kernels/INSADMomentumMeshAdvection.C
@@ -34,8 +34,8 @@ INSADMomentumMeshAdvection::validParams()
 
 INSADMomentumMeshAdvection::INSADMomentumMeshAdvection(const InputParameters & parameters)
   : ADVectorKernelValue(parameters),
-    _convected_mesh_strong_residual(
-        getADMaterialProperty<RealVectorValue>("convected_mesh_strong_residual"))
+    _advected_mesh_strong_residual(
+        getADMaterialProperty<RealVectorValue>("advected_mesh_strong_residual"))
 {
   setDisplacementParams(*this);
 }
@@ -43,5 +43,5 @@ INSADMomentumMeshAdvection::INSADMomentumMeshAdvection(const InputParameters & p
 ADRealVectorValue
 INSADMomentumMeshAdvection::precomputeQpResidual()
 {
-  return _convected_mesh_strong_residual[_qp];
+  return _advected_mesh_strong_residual[_qp];
 }
diff --git a/modules/navier_stokes/src/kernels/INSBase.C b/modules/navier_stokes/src/kernels/INSBase.C
index 4a182a542746..ffd096191a26 100644
--- a/modules/navier_stokes/src/kernels/INSBase.C
+++ b/modules/navier_stokes/src/kernels/INSBase.C
@@ -38,6 +38,9 @@ INSBase::validParams()
   params.addParam<bool>("transient_term",
                         false,
                         "Whether there should be a transient term in the momentum residuals.");
+  params.addCoupledVar("disp_x", "The x displacement");
+  params.addCoupledVar("disp_y", "The y displacement");
+  params.addCoupledVar("disp_z", "The z displacement");
 
   return params;
 }
@@ -88,10 +91,25 @@ INSBase::INSBase(const InputParameters & parameters)
     _alpha(getParam<Real>("alpha")),
     _laplace(getParam<bool>("laplace")),
     _convective_term(getParam<bool>("convective_term")),
-    _transient_term(getParam<bool>("transient_term"))
+    _transient_term(getParam<bool>("transient_term")),
+
+    // Displacements for mesh velocity for ALE simulations
+    _disps_provided(isParamValid("disp_x")),
+    _disp_x_dot(isParamValid("disp_x") ? coupledDot("disp_x") : _zero),
+    _disp_y_dot(isParamValid("disp_y") ? coupledDot("disp_y") : _zero),
+    _disp_z_dot(isParamValid("disp_z") ? coupledDot("disp_z") : _zero)
 {
 }
 
+RealVectorValue
+INSBase::relativeVelocity() const
+{
+  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
+  if (_disps_provided)
+    U -= RealVectorValue{_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]};
+  return U;
+}
+
 RealVectorValue
 INSBase::convectiveTerm()
 {
@@ -270,7 +288,7 @@ Real
 INSBase::tau()
 {
   Real nu = _mu[_qp] / _rho[_qp];
-  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
+  const auto U = relativeVelocity();
   Real h = _current_elem->hmax();
   Real transient_part = _transient_term ? 4. / (_dt * _dt) : 0.;
   return _alpha / std::sqrt(transient_part + (2. * U.norm() / h) * (2. * U.norm() / h) +
@@ -281,7 +299,7 @@ Real
 INSBase::tauNodal()
 {
   Real nu = _mu[_qp] / _rho[_qp];
-  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
+  const auto U = relativeVelocity();
   Real h = _current_elem->hmax();
   Real xi;
   if (nu < std::numeric_limits<Real>::epsilon())
@@ -298,7 +316,7 @@ Real
 INSBase::dTauDUComp(unsigned comp)
 {
   Real nu = _mu[_qp] / _rho[_qp];
-  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
+  const auto U = relativeVelocity();
   Real h = _current_elem->hmax();
   Real transient_part = _transient_term ? 4. / (_dt * _dt) : 0.;
   return -_alpha / 2. *
diff --git a/modules/navier_stokes/src/kernels/INSMomentumBase.C b/modules/navier_stokes/src/kernels/INSMomentumBase.C
index 59ea047a5794..35c605cdcc25 100644
--- a/modules/navier_stokes/src/kernels/INSMomentumBase.C
+++ b/modules/navier_stokes/src/kernels/INSMomentumBase.C
@@ -65,13 +65,11 @@ INSMomentumBase::computeQpResidual()
 Real
 INSMomentumBase::computeQpPGResidual()
 {
-  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
+  const auto U = relativeVelocity();
 
-  RealVectorValue convective_term = _convective_term ? convectiveTerm() : RealVectorValue(0, 0, 0);
-  RealVectorValue viscous_term =
-      _laplace ? strongViscousTermLaplace() : strongViscousTermTraction();
-  RealVectorValue transient_term =
-      _transient_term ? timeDerivativeTerm() : RealVectorValue(0, 0, 0);
+  const auto convective_term = _convective_term ? convectiveTerm() : RealVectorValue(0, 0, 0);
+  const auto viscous_term = _laplace ? strongViscousTermLaplace() : strongViscousTermTraction();
+  const auto transient_term = _transient_term ? timeDerivativeTerm() : RealVectorValue(0, 0, 0);
 
   return tau() * U * _grad_test[_i][_qp] *
          ((convective_term + viscous_term + transient_term + strongPressureTerm() +
@@ -109,18 +107,19 @@ INSMomentumBase::computeQpJacobian()
 Real
 INSMomentumBase::computeQpPGJacobian(unsigned comp)
 {
-  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
+  const auto U = relativeVelocity();
   RealVectorValue d_U_d_U_comp(0, 0, 0);
   d_U_d_U_comp(comp) = _phi[_j][_qp];
 
-  Real convective_term = _convective_term ? convectiveTerm()(_component) : 0;
-  Real d_convective_term_d_u_comp = _convective_term ? dConvecDUComp(comp)(_component) : 0;
-  Real viscous_term =
+  const Real convective_term = _convective_term ? convectiveTerm()(_component) : 0;
+  const Real d_convective_term_d_u_comp = _convective_term ? dConvecDUComp(comp)(_component) : 0;
+  const Real viscous_term =
       _laplace ? strongViscousTermLaplace()(_component) : strongViscousTermTraction()(_component);
-  Real d_viscous_term_d_u_comp = _laplace ? dStrongViscDUCompLaplace(comp)(_component)
-                                          : dStrongViscDUCompTraction(comp)(_component);
-  Real transient_term = _transient_term ? timeDerivativeTerm()(_component) : 0;
-  Real d_transient_term_d_u_comp = _transient_term ? dTimeDerivativeDUComp(comp)(_component) : 0;
+  const Real d_viscous_term_d_u_comp = _laplace ? dStrongViscDUCompLaplace(comp)(_component)
+                                                : dStrongViscDUCompTraction(comp)(_component);
+  const Real transient_term = _transient_term ? timeDerivativeTerm()(_component) : 0;
+  const Real d_transient_term_d_u_comp =
+      _transient_term ? dTimeDerivativeDUComp(comp)(_component) : 0;
 
   return dTauDUComp(comp) * U * _grad_test[_i][_qp] *
              (convective_term + viscous_term + strongPressureTerm()(_component) +
@@ -182,7 +181,7 @@ INSMomentumBase::computeQpOffDiagJacobian(unsigned jvar)
 
     if (_supg)
     {
-      RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
+      const auto U = relativeVelocity();
       jac += tau() * U * _grad_test[_i][_qp] * dStrongPressureDPressure()(_component);
     }
 
diff --git a/modules/navier_stokes/src/materials/INSAD3Eqn.C b/modules/navier_stokes/src/materials/INSAD3Eqn.C
index e64a5671f171..f5cd5d7f3c4a 100644
--- a/modules/navier_stokes/src/materials/INSAD3Eqn.C
+++ b/modules/navier_stokes/src/materials/INSAD3Eqn.C
@@ -36,8 +36,8 @@ INSAD3Eqn::INSAD3Eqn(const InputParameters & parameters)
         declareADProperty<Real>("temperature_ambient_convection_strong_residual")),
     _temperature_source_strong_residual(
         declareADProperty<Real>("temperature_source_strong_residual")),
-    _temperature_convected_mesh_strong_residual(
-        declareADProperty<Real>("temperature_convected_mesh_strong_residual"))
+    _temperature_advected_mesh_strong_residual(
+        declareADProperty<Real>("temperature_advected_mesh_strong_residual"))
 {
 }
 
@@ -127,7 +127,7 @@ INSAD3Eqn::computeQpProperties()
     }
   }
 
-  if (_has_convected_mesh)
-    _temperature_convected_mesh_strong_residual[_qp] =
+  if (_has_advected_mesh)
+    _temperature_advected_mesh_strong_residual[_qp] =
         -_rho[_qp] * _cp[_qp] * _mesh_velocity[_qp] * _grad_temperature[_qp];
 }
diff --git a/modules/navier_stokes/src/materials/INSADMaterial.C b/modules/navier_stokes/src/materials/INSADMaterial.C
index ae7b48313275..72a4291035c0 100644
--- a/modules/navier_stokes/src/materials/INSADMaterial.C
+++ b/modules/navier_stokes/src/materials/INSADMaterial.C
@@ -50,8 +50,9 @@ INSADMaterial::INSADMaterial(const InputParameters & parameters)
     _coupled_force_strong_residual(
         declareADProperty<RealVectorValue>("coupled_force_strong_residual")),
     _mesh_velocity(declareADProperty<RealVectorValue>("mesh_velocity")),
-    _convected_mesh_strong_residual(
-        declareADProperty<RealVectorValue>("convected_mesh_strong_residual")),
+    _relative_velocity(declareADProperty<RealVectorValue>("relative_velocity")),
+    _advected_mesh_strong_residual(
+        declareADProperty<RealVectorValue>("advected_mesh_strong_residual")),
     // _mms_function_strong_residual(declareProperty<RealVectorValue>("mms_function_strong_residual")),
     _use_displaced_mesh(getParam<bool>("use_displaced_mesh")),
     _ad_q_point(_bnd ? _assembly.adQPointsFace() : _assembly.adQPoints()),
@@ -112,8 +113,7 @@ INSADMaterial::subdomainSetup()
   // Setup data for Arbitrary Lagrangian Eulerian (ALE) simulations in which the simulation domain
   // is displacing. We will need to subtract the mesh velocity from the velocity solution in order
   // to get the correct material velocity for the momentum convection term.
-  if ((_has_convected_mesh =
-           _object_tracker->get<bool>("has_convected_mesh", _current_subdomain_id)))
+  if ((_has_advected_mesh = _object_tracker->get<bool>("has_advected_mesh", _current_subdomain_id)))
   {
     auto & disp_x = _subproblem.getStandardVariable(
         _tid, _object_tracker->get<VariableName>("disp_x", _current_subdomain_id));
@@ -210,16 +210,17 @@ INSADMaterial::computeQpProperties()
   if (_has_boussinesq)
     _boussinesq_strong_residual[_qp] = (*_boussinesq_alpha)[_qp] * _gravity_vector * _rho[_qp] *
                                        ((*_temperature)[_qp] - (*_ref_temp)[_qp]);
+  _relative_velocity[_qp] = _velocity[_qp];
 
-  if (_has_convected_mesh)
+  if (_has_advected_mesh)
   {
     _mesh_velocity[_qp](0) = (*_disp_x_dot)[_qp];
     if (_disp_y_dot)
       _mesh_velocity[_qp](1) = (*_disp_y_dot)[_qp];
     if (_disp_z_dot)
       _mesh_velocity[_qp](2) = (*_disp_z_dot)[_qp];
-    _convected_mesh_strong_residual[_qp] =
-        -_rho[_qp] * _grad_velocity[_qp].left_multiply(_mesh_velocity[_qp]);
+    _relative_velocity[_qp] -= _mesh_velocity[_qp];
+    _advected_mesh_strong_residual[_qp] = -_rho[_qp] * _grad_velocity[_qp] * _mesh_velocity[_qp];
   }
 
   if (_has_coupled_force)
diff --git a/modules/navier_stokes/src/materials/INSADStabilized3Eqn.C b/modules/navier_stokes/src/materials/INSADStabilized3Eqn.C
index e6690d7e3d9b..02e383d60082 100644
--- a/modules/navier_stokes/src/materials/INSADStabilized3Eqn.C
+++ b/modules/navier_stokes/src/materials/INSADStabilized3Eqn.C
@@ -46,9 +46,8 @@ INSADStabilized3Eqn::computeQpProperties()
 
   const auto dissipation_coefficient = _k[_qp] / (_rho[_qp] * _cp[_qp]);
   const auto transient_part = _has_energy_transient ? 4. / (_dt * _dt) : 0.;
-  const auto speed = NS::computeSpeed(_velocity[_qp]);
   _tau_energy[_qp] =
-      _alpha / std::sqrt(transient_part + (2. * speed / _hmax) * (2. * speed / _hmax) +
+      _alpha / std::sqrt(transient_part + (2. * _speed / _hmax) * (2. * _speed / _hmax) +
                          9. * (4. * dissipation_coefficient / (_hmax * _hmax)) *
                              (4. * dissipation_coefficient / (_hmax * _hmax)));
 
@@ -69,6 +68,6 @@ INSADStabilized3Eqn::computeQpProperties()
   if (_has_heat_source)
     _temperature_strong_residual[_qp] += _temperature_source_strong_residual[_qp];
 
-  if (_has_convected_mesh)
-    _temperature_strong_residual[_qp] += _temperature_convected_mesh_strong_residual[_qp];
+  if (_has_advected_mesh)
+    _temperature_strong_residual[_qp] += _temperature_advected_mesh_strong_residual[_qp];
 }
diff --git a/modules/navier_stokes/src/userobjects/INSADObjectTracker.C b/modules/navier_stokes/src/userobjects/INSADObjectTracker.C
index 0eef5ca4a39f..746cab509dbf 100644
--- a/modules/navier_stokes/src/userobjects/INSADObjectTracker.C
+++ b/modules/navier_stokes/src/userobjects/INSADObjectTracker.C
@@ -62,14 +62,14 @@ INSADObjectTracker::validTrackerParams()
   params.addParam<std::vector<FunctionName>>("coupled_force_vector_function",
                                              "The function(s) standing in as a coupled force(s)");
 
-  params.addParam<bool>("has_convected_mesh",
+  params.addParam<bool>("has_advected_mesh",
                         false,
                         "Whether the fluid domain is undergoing displacement in which case we must "
-                        "add a convecting mesh term to correct the material velocity.");
+                        "add an advecting mesh term to correct the material velocity.");
   params.addParam<VariableName>("disp_x", "The x displacement");
   params.addParam<VariableName>("disp_y", "The y displacement");
   params.addParam<VariableName>("disp_z", "The z displacement");
-  params.addParamNamesToGroup("has_convected_mesh disp_x disp_y disp_z", "Moving mesh");
+  params.addParamNamesToGroup("has_advected_mesh disp_x disp_y disp_z", "Moving mesh");
   return params;
 }