diff --git a/modules/navier_stokes/include/kernels/PINSFEFluidTemperatureTimeDerivative.h b/modules/navier_stokes/include/kernels/PINSFEFluidTemperatureTimeDerivative.h
index 9f220d8c73a9..05f6a26de574 100644
--- a/modules/navier_stokes/include/kernels/PINSFEFluidTemperatureTimeDerivative.h
+++ b/modules/navier_stokes/include/kernels/PINSFEFluidTemperatureTimeDerivative.h
@@ -27,10 +27,13 @@ class PINSFEFluidTemperatureTimeDerivative : public TimeDerivative
 protected:
   virtual Real computeQpResidual() override;
   virtual Real computeQpJacobian() override;
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;
 
   bool _conservative_form;
   const VariableValue & _pressure;
   const VariableValue & _pressure_dot;
+  const VariableValue & _d_pressure_dot_du;
+  unsigned int _pressure_var_number;
   const VariableValue & _porosity;
   const MaterialProperty<Real> & _rho;
   const MaterialProperty<Real> & _cp;
diff --git a/modules/navier_stokes/include/kernels/PINSFEFluidVelocityTimeDerivative.h b/modules/navier_stokes/include/kernels/PINSFEFluidVelocityTimeDerivative.h
index 7b4a2f284076..32a259394167 100644
--- a/modules/navier_stokes/include/kernels/PINSFEFluidVelocityTimeDerivative.h
+++ b/modules/navier_stokes/include/kernels/PINSFEFluidVelocityTimeDerivative.h
@@ -25,14 +25,19 @@ class PINSFEFluidVelocityTimeDerivative : public TimeDerivative
   PINSFEFluidVelocityTimeDerivative(const InputParameters & parameters);
 
 protected:
-  virtual Real computeQpResidual();
-  virtual Real computeQpJacobian();
+  virtual Real computeQpResidual() override;
+  virtual Real computeQpJacobian() override;
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;
 
   bool _conservative_form;
   const VariableValue & _pressure;
   const VariableValue & _temperature;
   const VariableValue & _temperature_dot;
   const VariableValue & _pressure_dot;
+  const VariableValue & _d_pressure_dot_du;
+  const VariableValue & _d_temperature_dot_du;
+  unsigned int _pressure_var_number;
+  unsigned int _temperature_var_number;
   const MaterialProperty<Real> & _rho;
   const SinglePhaseFluidProperties & _eos;
 };
diff --git a/modules/navier_stokes/src/kernels/INSFEFluidEnergyKernel.C b/modules/navier_stokes/src/kernels/INSFEFluidEnergyKernel.C
index cd347c9c4def..997e6a425f1e 100644
--- a/modules/navier_stokes/src/kernels/INSFEFluidEnergyKernel.C
+++ b/modules/navier_stokes/src/kernels/INSFEFluidEnergyKernel.C
@@ -138,6 +138,28 @@ INSFEFluidEnergyKernel::computeQpOffDiagJacobian(unsigned int jvar)
   Real jac = 0.;
   switch (m)
   {
+    case 0:
+    {
+      // pressure variable: differentiate convection and SUPG terms through EOS (rho, h)
+      Real porosity = _has_porosity ? _porosity[_qp] : 1.0;
+      Real rho, drho_dp, drho_dT;
+      _eos.rho_from_p_T(_pressure[_qp], _u[_qp], rho, drho_dp, drho_dT);
+      if (_conservative_form)
+      {
+        Real h, dh_dp, dh_dT;
+        _eos.h_from_p_T(_pressure[_qp], _u[_qp], h, dh_dp, dh_dT);
+        jac -= (drho_dp * h + rho * dh_dp) * _phi[_j][_qp] * vec_vel * _grad_test[_i][_qp];
+      }
+      else
+        jac += drho_dp * _cp[_qp] * _phi[_j][_qp] * vec_vel * _grad_u[_qp] * _test[_i][_qp];
+      // SUPG contribution (strong form uses non-conservative convection regardless of form)
+      Real transient_supg = _bTransient ? porosity * drho_dp * _cp[_qp] * _u_dot[_qp] : 0.0;
+      Real convection_supg = drho_dp * _cp[_qp] * vec_vel * _grad_u[_qp];
+      jac += _phi[_j][_qp] * _taue[_qp] * _vel_elem * _grad_test[_i][_qp] *
+             (transient_supg + convection_supg);
+      break;
+    }
+
     case 1:
     case 2:
     case 3:
diff --git a/modules/navier_stokes/src/kernels/PINSFEFluidTemperatureTimeDerivative.C b/modules/navier_stokes/src/kernels/PINSFEFluidTemperatureTimeDerivative.C
index 6128191123ba..257815b84acb 100644
--- a/modules/navier_stokes/src/kernels/PINSFEFluidTemperatureTimeDerivative.C
+++ b/modules/navier_stokes/src/kernels/PINSFEFluidTemperatureTimeDerivative.C
@@ -35,6 +35,8 @@ PINSFEFluidTemperatureTimeDerivative::PINSFEFluidTemperatureTimeDerivative(
     _conservative_form(getParam<bool>("conservative_form")),
     _pressure(coupledValue("pressure")),
     _pressure_dot(coupledDot("pressure")),
+    _d_pressure_dot_du(coupledDotDu("pressure")),
+    _pressure_var_number(coupled("pressure")),
     _porosity(coupledValue("porosity")),
     _rho(getMaterialProperty<Real>("rho_fluid")),
     _cp(getMaterialProperty<Real>("cp_fluid")),
@@ -66,8 +68,26 @@ PINSFEFluidTemperatureTimeDerivative::computeQpJacobian()
     Real rho, drho_dp, drho_dT;
     _eos.rho_from_p_T(_pressure[_qp], _u[_qp], rho, drho_dp, drho_dT);
     Real drho_dt = drho_dT * _u_dot[_qp] + drho_dp * _pressure_dot[_qp];
+    // d(u)/dT = phi_j contribution
     jac += _porosity[_qp] * _cp[_qp] * _phi[_j][_qp] * drho_dt * _test[_i][_qp];
+    // d(drho_dt)/dT = drho_dT * d(T_dot)/dT_j contribution
+    jac += _porosity[_qp] * _cp[_qp] * _u[_qp] * drho_dT * _du_dot_du[_qp] * _phi[_j][_qp] *
+           _test[_i][_qp];
   }
 
   return jac;
 }
+
+Real
+PINSFEFluidTemperatureTimeDerivative::computeQpOffDiagJacobian(unsigned int jvar)
+{
+  if (_conservative_form && jvar == _pressure_var_number)
+  {
+    Real rho, drho_dp, drho_dT;
+    _eos.rho_from_p_T(_pressure[_qp], _u[_qp], rho, drho_dp, drho_dT);
+    // d(drho_dt)/dp = drho_dp * d(p_dot)/dp_j contribution from the conservative term
+    return _porosity[_qp] * _cp[_qp] * _u[_qp] * drho_dp * _d_pressure_dot_du[_qp] * _phi[_j][_qp] *
+           _test[_i][_qp];
+  }
+  return 0.0;
+}
diff --git a/modules/navier_stokes/src/kernels/PINSFEFluidVelocityTimeDerivative.C b/modules/navier_stokes/src/kernels/PINSFEFluidVelocityTimeDerivative.C
index b0af358f81ec..7783e0ad04b4 100644
--- a/modules/navier_stokes/src/kernels/PINSFEFluidVelocityTimeDerivative.C
+++ b/modules/navier_stokes/src/kernels/PINSFEFluidVelocityTimeDerivative.C
@@ -37,6 +37,10 @@ PINSFEFluidVelocityTimeDerivative::PINSFEFluidVelocityTimeDerivative(
     _temperature(coupledValue("temperature")),
     _temperature_dot(coupledDot("temperature")),
     _pressure_dot(coupledDot("pressure")),
+    _d_pressure_dot_du(coupledDotDu("pressure")),
+    _d_temperature_dot_du(coupledDotDu("temperature")),
+    _pressure_var_number(coupled("pressure")),
+    _temperature_var_number(coupled("temperature")),
     _rho(getMaterialProperty<Real>("rho_fluid")),
     _eos(getUserObject<SinglePhaseFluidProperties>("eos"))
 {
@@ -71,3 +75,18 @@ PINSFEFluidVelocityTimeDerivative::computeQpJacobian()
 
   return jac;
 }
+
+Real
+PINSFEFluidVelocityTimeDerivative::computeQpOffDiagJacobian(unsigned int jvar)
+{
+  if (_conservative_form)
+  {
+    Real rho, drho_dp, drho_dT;
+    _eos.rho_from_p_T(_pressure[_qp], _temperature[_qp], rho, drho_dp, drho_dT);
+    if (jvar == _pressure_var_number)
+      return _u[_qp] * drho_dp * _d_pressure_dot_du[_qp] * _phi[_j][_qp] * _test[_i][_qp];
+    else if (jvar == _temperature_var_number)
+      return _u[_qp] * drho_dT * _d_temperature_dot_du[_qp] * _phi[_j][_qp] * _test[_i][_qp];
+  }
+  return 0.0;
+}
diff --git a/modules/navier_stokes/src/materials/INSFEMaterial.C b/modules/navier_stokes/src/materials/INSFEMaterial.C
index cf2797c0460a..5067cc555fc6 100644
--- a/modules/navier_stokes/src/materials/INSFEMaterial.C
+++ b/modules/navier_stokes/src/materials/INSFEMaterial.C
@@ -212,4 +212,12 @@ INSFEMaterial::computeTau()
   _tauc[_qp] = 1. / std::sqrt(sqrt_term_pspg);
   _taum[_qp] = 1. / std::sqrt(sqrt_term_supg + visc_term * visc_term);
   _taue[_qp] = 1. / std::sqrt(sqrt_term_supg + diff_term * diff_term);
+
+  // add lower and upper bounds
+  Real tau_base = std::min(_hsupg[_qp] / (2 * _scaling_velocity), _dt / 2);
+  _tauc[_qp] = std::clamp(_tauc[_qp], 0.005 * tau_base, tau_base);
+
+  tau_base = std::min(_hsupg[_qp] / (2 * _vel_mag), _dt / 2);
+  _taum[_qp] = std::clamp(_taum[_qp], 0.01 * tau_base, 2.0 * tau_base);
+  _taue[_qp] = std::clamp(_taue[_qp], 0.01 * tau_base, 2.0 * tau_base);
 }
diff --git a/modules/navier_stokes/test/tests/finite_element/ins/lid_driven/gold/lid_driven_md_out.e b/modules/navier_stokes/test/tests/finite_element/ins/lid_driven/gold/lid_driven_md_out.e
index 3302fc6ce35e..f31e341dad10 100644
Binary files a/modules/navier_stokes/test/tests/finite_element/ins/lid_driven/gold/lid_driven_md_out.e and b/modules/navier_stokes/test/tests/finite_element/ins/lid_driven/gold/lid_driven_md_out.e differ
diff --git a/modules/navier_stokes/test/tests/finite_element/pins/expansion-channel/gold/expansion-channel-slip-wall_out.e b/modules/navier_stokes/test/tests/finite_element/pins/expansion-channel/gold/expansion-channel-slip-wall_out.e
index 6731b42a15ce..607418443daa 100644
Binary files a/modules/navier_stokes/test/tests/finite_element/pins/expansion-channel/gold/expansion-channel-slip-wall_out.e and b/modules/navier_stokes/test/tests/finite_element/pins/expansion-channel/gold/expansion-channel-slip-wall_out.e differ