diff --git a/framework/include/systems/DisplacedSystem.h b/framework/include/systems/DisplacedSystem.h
index 436e97404a28..00ffcc6a5c0c 100644
--- a/framework/include/systems/DisplacedSystem.h
+++ b/framework/include/systems/DisplacedSystem.h
@@ -149,9 +149,9 @@ class DisplacedSystem : public SystemBase
     return _undisplaced_system.solutionUDotDotOld();
   }
 
-  virtual Number & duDotDu() override { return _undisplaced_system.duDotDu(); }
+  virtual std::vector<Number> & duDotDu() override { return _undisplaced_system.duDotDu(); }
   virtual Number & duDotDotDu() override { return _undisplaced_system.duDotDotDu(); }
-  virtual const Number & duDotDu() const override { return _undisplaced_system.duDotDu(); }
+  virtual const Number & duDotDu(const unsigned int var_num) const override;
   virtual const Number & duDotDotDu() const override { return _undisplaced_system.duDotDotDu(); }
 
   virtual void addDotVectors() override { _undisplaced_system.addDotVectors(); }
@@ -294,3 +294,9 @@ DisplacedSystem::hasSolutionState(const unsigned int state,
 {
   return _undisplaced_system.hasSolutionState(state, iteration_type);
 }
+
+inline const Number &
+DisplacedSystem::duDotDu(const unsigned int var_num) const
+{
+  return _undisplaced_system.duDotDu(var_num);
+}
diff --git a/framework/include/systems/SystemBase.h b/framework/include/systems/SystemBase.h
index 1012e3e7dbd3..7d9e6e8b174e 100644
--- a/framework/include/systems/SystemBase.h
+++ b/framework/include/systems/SystemBase.h
@@ -231,9 +231,9 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
    */
   virtual void addDotVectors();
 
-  virtual Number & duDotDu() { return _du_dot_du; }
+  virtual std::vector<Number> & duDotDu() { return _du_dot_du; }
   virtual Number & duDotDotDu() { return _du_dotdot_du; }
-  virtual const Number & duDotDu() const { return _du_dot_du; }
+  virtual const Number & duDotDu(const unsigned int var_num) const { return _du_dot_du[var_num]; }
   virtual const Number & duDotDotDu() const { return _du_dotdot_du; }
 
   virtual NumericVector<Number> * solutionUDot() { return _u_dot; }
@@ -995,7 +995,9 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
   /// old solution vector for u^dotdot
   NumericVector<Number> * _u_dotdot_old;
 
-  Real _du_dot_du;
+  /// Derivative of time derivative of u with respect to uj. This depends on the time integration
+  /// scheme
+  std::vector<Real> _du_dot_du;
   Real _du_dotdot_du;
 
   /// Tagged vectors (pointer)
diff --git a/framework/include/timeintegrators/TimeIntegrator.h b/framework/include/timeintegrators/TimeIntegrator.h
index 99bb604ac1d3..b6d7bfd3ab9a 100644
--- a/framework/include/timeintegrators/TimeIntegrator.h
+++ b/framework/include/timeintegrators/TimeIntegrator.h
@@ -198,11 +198,15 @@ class TimeIntegrator : public MooseObject, public Restartable
   /// residual vector for non-time contributions
   NumericVector<Number> & _Re_non_time;
 
-  /// Derivative of time derivative with respect to current solution: \f$ {du^dot}\over{du} \f$
-  Real & _du_dot_du;
+  /// Derivative of time derivative with respect to current solution: \f$ {du^dot}\over{du} \f$ for
+  /// the different variables. We will only modify the elements in this vector corresponding to the
+  /// variables that we integrate
+  std::vector<Real> & _du_dot_du;
   /// solution vectors
   const NumericVector<Number> * const & _solution;
   const NumericVector<Number> & _solution_old;
+  std::unique_ptr<NumericVector<Number>> _solution_sub;
+  std::unique_ptr<NumericVector<Number>> _solution_old_sub;
   //
   int & _t_step;
   //
diff --git a/framework/src/systems/SystemBase.C b/framework/src/systems/SystemBase.C
index afba51b495bd..6f379af3301f 100644
--- a/framework/src/systems/SystemBase.C
+++ b/framework/src/systems/SystemBase.C
@@ -800,6 +800,7 @@ SystemBase::addVariable(const std::string & var_type,
   // getMaxVariableNumber is an API method used in Rattlesnake
   if (var_num > _max_var_number)
     _max_var_number = var_num;
+  _du_dot_du.resize(var_num + 1);
 }
 
 bool
diff --git a/framework/src/timeintegrators/AStableDirk4.C b/framework/src/timeintegrators/AStableDirk4.C
index 5e6cd0af00cd..4a23f5580092 100644
--- a/framework/src/timeintegrators/AStableDirk4.C
+++ b/framework/src/timeintegrators/AStableDirk4.C
@@ -92,7 +92,9 @@ AStableDirk4::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-  _du_dot_du = 1. / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1. / _dt;
 }
 
 void
diff --git a/framework/src/timeintegrators/ActuallyExplicitEuler.C b/framework/src/timeintegrators/ActuallyExplicitEuler.C
index 0d3f332a16e4..5b4a154da796 100644
--- a/framework/src/timeintegrators/ActuallyExplicitEuler.C
+++ b/framework/src/timeintegrators/ActuallyExplicitEuler.C
@@ -50,7 +50,9 @@ ActuallyExplicitEuler::computeTimeDerivatives()
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
 
-  _du_dot_du = 1.0 / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1.0 / _dt;
 }
 
 void
diff --git a/framework/src/timeintegrators/BDF2.C b/framework/src/timeintegrators/BDF2.C
index 416961d7f432..43a3f6bb15ab 100644
--- a/framework/src/timeintegrators/BDF2.C
+++ b/framework/src/timeintegrators/BDF2.C
@@ -52,12 +52,16 @@ BDF2::computeTimeDerivatives()
   if (_t_step == 1)
   {
     u_dot = *_solution;
-    _du_dot_du = 1. / _dt;
+    for (const auto i : index_range(_du_dot_du))
+      if (integratesVar(i))
+        _du_dot_du[i] = 1. / _dt;
   }
   else
   {
     u_dot.zero();
-    _du_dot_du = _weight[0] / _dt;
+    for (const auto i : index_range(_du_dot_du))
+      if (integratesVar(i))
+        _du_dot_du[i] = _weight[0] / _dt;
   }
   computeTimeDerivativeHelper(u_dot, *_solution, _solution_old, _solution_older);
   u_dot.close();
diff --git a/framework/src/timeintegrators/CentralDifference.C b/framework/src/timeintegrators/CentralDifference.C
index a3cfb30eff89..869df4988321 100644
--- a/framework/src/timeintegrators/CentralDifference.C
+++ b/framework/src/timeintegrators/CentralDifference.C
@@ -87,7 +87,9 @@ CentralDifference::computeTimeDerivatives()
   u_dot.close();
 
   // used for Jacobian calculations
-  _du_dot_du = 1.0 / (2 * _dt);
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1.0 / (2 * _dt);
   _du_dotdot_du = 1.0 / (_dt * _dt);
 
   // Computing udotdot "factor"
diff --git a/framework/src/timeintegrators/CrankNicolson.C b/framework/src/timeintegrators/CrankNicolson.C
index be4f3c7ba631..4fcf5d149bbb 100644
--- a/framework/src/timeintegrators/CrankNicolson.C
+++ b/framework/src/timeintegrators/CrankNicolson.C
@@ -34,11 +34,26 @@ CrankNicolson::computeTimeDerivatives()
                "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
 
   NumericVector<Number> & u_dot = *_sys.solutionUDot();
-  u_dot = *_solution;
-  computeTimeDerivativeHelper(u_dot, _solution_old);
-  u_dot.close();
+  if (!_var_restriction)
+  {
+    u_dot = *_solution;
+    computeTimeDerivativeHelper(u_dot, _solution_old);
+  }
+  else
+  {
+    auto u_dot_sub = u_dot.get_subvector(_local_indices);
+    _solution->create_subvector(*_solution_sub, _local_indices, false);
+    _solution_old.create_subvector(*_solution_old_sub, _local_indices, false);
+    *u_dot_sub = *_solution_sub;
+    computeTimeDerivativeHelper(*u_dot_sub, *_solution_old_sub);
+    u_dot.restore_subvector(std::move(*u_dot_sub), _local_indices);
+    // Scatter info needed for ghosts
+    u_dot.close();
+  }
 
-  _du_dot_du = 2. / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 2. / _dt;
 }
 
 void
@@ -61,7 +76,9 @@ CrankNicolson::init()
   // time derivative is assumed to be zero on initial
   NumericVector<Number> & u_dot = *_sys.solutionUDot();
   u_dot.zero();
-  _du_dot_du = 0;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 0;
 
   // compute residual for the initial time step
   // Note: we can not directly pass _residual_old in computeResidualTag because
diff --git a/framework/src/timeintegrators/ExplicitEuler.C b/framework/src/timeintegrators/ExplicitEuler.C
index 1639be7cc1ad..bf56516eae59 100644
--- a/framework/src/timeintegrators/ExplicitEuler.C
+++ b/framework/src/timeintegrators/ExplicitEuler.C
@@ -44,7 +44,9 @@ ExplicitEuler::computeTimeDerivatives()
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
 
-  _du_dot_du = 1.0 / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1.0 / _dt;
 }
 
 void
diff --git a/framework/src/timeintegrators/ExplicitRK2.C b/framework/src/timeintegrators/ExplicitRK2.C
index 8f8d95889463..b99af6045e1c 100644
--- a/framework/src/timeintegrators/ExplicitRK2.C
+++ b/framework/src/timeintegrators/ExplicitRK2.C
@@ -54,7 +54,9 @@ ExplicitRK2::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old, _solution_older);
 
-  _du_dot_du = 1. / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1. / _dt;
   u_dot.close();
 }
 
diff --git a/framework/src/timeintegrators/ExplicitSSPRungeKutta.C b/framework/src/timeintegrators/ExplicitSSPRungeKutta.C
index 9e1b54196188..f238db646c4b 100644
--- a/framework/src/timeintegrators/ExplicitSSPRungeKutta.C
+++ b/framework/src/timeintegrators/ExplicitSSPRungeKutta.C
@@ -69,7 +69,9 @@ void
 ExplicitSSPRungeKutta::computeTimeDerivatives()
 {
   // Only the Jacobian needs to be computed, since the mass matrix needs it
-  _du_dot_du = 1.0 / (_b[_stage] * _dt);
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1.0 / (_b[_stage] * _dt);
 }
 
 void
diff --git a/framework/src/timeintegrators/ExplicitTVDRK2.C b/framework/src/timeintegrators/ExplicitTVDRK2.C
index 1a625c6d904b..aed081626516 100644
--- a/framework/src/timeintegrators/ExplicitTVDRK2.C
+++ b/framework/src/timeintegrators/ExplicitTVDRK2.C
@@ -56,7 +56,9 @@ ExplicitTVDRK2::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old, _solution_older);
 
-  _du_dot_du = 1. / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1. / _dt;
   u_dot.close();
 }
 
diff --git a/framework/src/timeintegrators/ImplicitEuler.C b/framework/src/timeintegrators/ImplicitEuler.C
index ba3aba153729..7cf2c2af486a 100644
--- a/framework/src/timeintegrators/ImplicitEuler.C
+++ b/framework/src/timeintegrators/ImplicitEuler.C
@@ -32,11 +32,26 @@ ImplicitEuler::computeTimeDerivatives()
                "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
 
   NumericVector<Number> & u_dot = *_sys.solutionUDot();
-  u_dot = *_solution;
-  computeTimeDerivativeHelper(u_dot, _solution_old);
-  u_dot.close();
+  if (!_var_restriction)
+  {
+    u_dot = *_solution;
+    computeTimeDerivativeHelper(u_dot, _solution_old);
+  }
+  else
+  {
+    auto u_dot_sub = u_dot.get_subvector(_local_indices);
+    _solution->create_subvector(*_solution_sub, _local_indices, false);
+    _solution_old.create_subvector(*_solution_old_sub, _local_indices, false);
+    *u_dot_sub = *_solution_sub;
+    computeTimeDerivativeHelper(*u_dot_sub, *_solution_old_sub);
+    u_dot.restore_subvector(std::move(*u_dot_sub), _local_indices);
+    // Scatter info needed for ghosts
+    u_dot.close();
+  }
 
-  _du_dot_du = 1.0 / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1.0 / _dt;
 }
 
 void
diff --git a/framework/src/timeintegrators/ImplicitMidpoint.C b/framework/src/timeintegrators/ImplicitMidpoint.C
index 7d7daa654ebe..8bb834dec8de 100644
--- a/framework/src/timeintegrators/ImplicitMidpoint.C
+++ b/framework/src/timeintegrators/ImplicitMidpoint.C
@@ -45,7 +45,9 @@ ImplicitMidpoint::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-  _du_dot_du = 1. / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1. / _dt;
 }
 
 void
diff --git a/framework/src/timeintegrators/LStableDirk2.C b/framework/src/timeintegrators/LStableDirk2.C
index 5904203e5ea2..cfd22dc8d25b 100644
--- a/framework/src/timeintegrators/LStableDirk2.C
+++ b/framework/src/timeintegrators/LStableDirk2.C
@@ -48,7 +48,9 @@ LStableDirk2::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-  _du_dot_du = 1. / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1. / _dt;
 }
 
 void
diff --git a/framework/src/timeintegrators/LStableDirk3.C b/framework/src/timeintegrators/LStableDirk3.C
index 9fec25a0deca..0412e0ae680b 100644
--- a/framework/src/timeintegrators/LStableDirk3.C
+++ b/framework/src/timeintegrators/LStableDirk3.C
@@ -67,7 +67,9 @@ LStableDirk3::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-  _du_dot_du = 1. / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1. / _dt;
 }
 
 void
diff --git a/framework/src/timeintegrators/LStableDirk4.C b/framework/src/timeintegrators/LStableDirk4.C
index f144ced1ee8b..ea9af65567da 100644
--- a/framework/src/timeintegrators/LStableDirk4.C
+++ b/framework/src/timeintegrators/LStableDirk4.C
@@ -62,7 +62,9 @@ LStableDirk4::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-  _du_dot_du = 1. / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = 1. / _dt;
 }
 
 void
diff --git a/framework/src/timeintegrators/NewmarkBeta.C b/framework/src/timeintegrators/NewmarkBeta.C
index 4fa100cbdef1..0fc03741a4f5 100644
--- a/framework/src/timeintegrators/NewmarkBeta.C
+++ b/framework/src/timeintegrators/NewmarkBeta.C
@@ -92,7 +92,9 @@ NewmarkBeta::computeTimeDerivatives()
 
   // used for Jacobian calculations
   _du_dotdot_du = 1.0 / _beta / _dt / _dt;
-  _du_dot_du = _gamma / _beta / _dt;
+  for (const auto i : index_range(_du_dot_du))
+    if (integratesVar(i))
+      _du_dot_du[i] = _gamma / _beta / _dt;
 }
 
 void
diff --git a/framework/src/timeintegrators/TimeIntegrator.C b/framework/src/timeintegrators/TimeIntegrator.C
index 311e94eca937..b8451b816e64 100644
--- a/framework/src/timeintegrators/TimeIntegrator.C
+++ b/framework/src/timeintegrators/TimeIntegrator.C
@@ -89,6 +89,9 @@ TimeIntegrator::init()
                var_dof_indices.end(),
                std::back_inserter(_local_indices));
   }
+
+  _solution_sub = NumericVector<Number>::build(_solution->comm());
+  _solution_old_sub = NumericVector<Number>::build(_solution_old.comm());
 }
 
 void
diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C
index 3458868d59e9..39e179fc9ecd 100644
--- a/framework/src/variables/MooseVariableData.C
+++ b/framework/src/variables/MooseVariableData.C
@@ -1027,7 +1027,7 @@ MooseVariableData<OutputType>::computeMonomialValues()
   Real u_dotdot = 0;
   Real u_dot_old = 0;
   Real u_dotdot_old = 0;
-  const Real & du_dot_du = _sys.duDotDu();
+  const Real & du_dot_du = _sys.duDotDu(_var_num);
   const Real & du_dotdot_du = _sys.duDotDotDu();
 
   if (is_transient)
diff --git a/framework/src/variables/MooseVariableDataBase.C b/framework/src/variables/MooseVariableDataBase.C
index 34917c824b06..da3c3fc96029 100644
--- a/framework/src/variables/MooseVariableDataBase.C
+++ b/framework/src/variables/MooseVariableDataBase.C
@@ -600,7 +600,7 @@ MooseVariableDataBase<OutputType>::fetchDoFValues()
   {
     _dof_du_dot_du.resize(n);
     for (decltype(n) i = 0; i < n; ++i)
-      _dof_du_dot_du[i] = _sys.duDotDu();
+      _dof_du_dot_du[i] = _sys.duDotDu(_var.number());
   }
   if (_need_du_dotdot_du || _need_dof_du_dotdot_du)
   {
@@ -708,7 +708,7 @@ MooseVariableDataBase<RealEigenVector>::fetchDoFValues()
   {
     _dof_du_dot_du.resize(n);
     for (decltype(n) i = 0; i < n; ++i)
-      _dof_du_dot_du[i] = _sys.duDotDu();
+      _dof_du_dot_du[i] = _sys.duDotDu(_var.number());
   }
   if (_need_du_dotdot_du || _need_dof_du_dotdot_du)
   {
diff --git a/framework/src/variables/MooseVariableScalar.C b/framework/src/variables/MooseVariableScalar.C
index 9779c8816f54..f8a9616a9916 100644
--- a/framework/src/variables/MooseVariableScalar.C
+++ b/framework/src/variables/MooseVariableScalar.C
@@ -100,7 +100,7 @@ MooseVariableScalar::reinit(bool reinit_for_derivative_reordering /* = false*/)
   const NumericVector<Real> * u_dotdot = sys.solutionUDotDot();
   const NumericVector<Real> * u_dot_old = sys.solutionUDotOld();
   const NumericVector<Real> * u_dotdot_old = sys.solutionUDotDotOld();
-  const Real & du_dot_du = sys.duDotDu();
+  const Real & du_dot_du = sys.duDotDu(number());
   const Real & du_dotdot_du = sys.duDotDotDu();
   auto safe_access_tagged_vectors = sys.subproblem().safeAccessTaggedVectors();
   auto safe_access_tagged_matrices = sys.subproblem().safeAccessTaggedMatrices();