diff --git a/framework/doc/content/source/linearfvkernels/LinearFVTimeDerivative.md b/framework/doc/content/source/linearfvkernels/LinearFVTimeDerivative.md
new file mode 100644
index 000000000000..62482e14afb0
--- /dev/null
+++ b/framework/doc/content/source/linearfvkernels/LinearFVTimeDerivative.md
@@ -0,0 +1,47 @@
+# LinearFVTimeDerivative
+
+## Description
+
+This kernel represents a time derivative term in a partial differential equation
+discretized using the finite volume method:
+
+!equation
+\int\limits_{V_C} \frac{\partial cu}{\partial t} dV \approx \left(\frac{\partial cu}{\partial t}\right)_C V_C~,
+
+where $\left(\frac{\partial cu}{\partial t}\right)_C$ and $V_C$ are the time derivative of the
+field value at the cell center and the cell volume, respectively.
+Note that we added a multiplier, $c$ which often represents a material property.
+A good example for the multiplier can be the density in the momentum equation
+in the Navier Stokes equation.
+This can be defined through parameter [!param](/LinearFVKernels/LinearFVTimeDerivative/factor)
+that accepts anything that supports functor-based evaluations. For more information on functors in
+MOOSE, see [Functors/index.md].
+This kernel adds to the matrix diagonal and right hand side of a
+linear system and the contributions depend on the
+method chosen for time integration. For more information on available methods, see
+the [TimeIntegrators](Executioner/TimeIntegrators/index.md) page.
+For example, with an implicit Euler scheme the contribution to the right hand side becomes:
+
+!equation
+\frac{c_C}{\Delta t}V_C,
+
+where $\Delta t$ and $c_C$ are the time step size and multiplier at the cell center,
+respectively. With these, the contribution to the right hand side becomes:
+
+!equation
+\frac{c_C u_{old,C}}{\Delta t}V_C,
+
+where $u_{old,C}$ represents the solution at the previous time step.
+
+## Example Syntax
+
+The case below demonstrates the use of `LinearFVTimeDerivative` used in a simple
+linear time-dependent diffusion problem:
+
+!listing test/tests/time_integrators/implicit-euler/ie-linearfv.i block=LinearFVKernels
+
+!syntax parameters /LinearFVKernels/LinearFVTimeDerivative
+
+!syntax inputs /LinearFVKernels/LinearFVTimeDerivative
+
+!syntax children /LinearFVKernels/LinearFVTimeDerivative
diff --git a/framework/doc/content/source/timeintegrators/ImplicitEuler.md b/framework/doc/content/source/timeintegrators/ImplicitEuler.md
index 6648a2cfeb90..83a4378ede34 100644
--- a/framework/doc/content/source/timeintegrators/ImplicitEuler.md
+++ b/framework/doc/content/source/timeintegrators/ImplicitEuler.md
@@ -21,6 +21,22 @@ This is an implicit system with $U(t+\Delta t)$, the variable to solve for, appe
 equation. We solve this system iteratively, usually with a Newton or Newton-Krylov method as described
 in the non linear system solve [documentation](systems/NonlinearSystem.md).
 
+## Contributions to linear systems
+
+For [linear systems](systems/LinearSystem.md), on top of creating the
+time derivatives of the degrees of freedom, this provides contributions
+to the matrix diagonal and the right hand side. Taking a finite volume system for example,
+the contributions to the matrix diagonal will be:
+
+!equation
+\frac{1}{\Delta t}V_C,
+
+where $\Delta t$ and $\frac{V_C}$ are the time step size and cell volume,
+respectively. The contribution to the right hand side is:
+
+!equation
+\frac{u_{old,C}}{\Delta t}V_C,
+
 !syntax parameters /Executioner/TimeIntegrator/ImplicitEuler
 
 !syntax inputs /Executioner/TimeIntegrator/ImplicitEuler
diff --git a/framework/include/linearfvkernels/LinearFVElementalKernel.h b/framework/include/linearfvkernels/LinearFVElementalKernel.h
index b9af62b54dd8..0c1e08f490df 100644
--- a/framework/include/linearfvkernels/LinearFVElementalKernel.h
+++ b/framework/include/linearfvkernels/LinearFVElementalKernel.h
@@ -35,10 +35,11 @@ class LinearFVElementalKernel : public LinearFVKernel
    * Set the current ElemInfo object.
    * @param elem_info The new ElemInfo which will be used as the current
    */
-  void setCurrentElemInfo(const ElemInfo * elem_info) { _current_elem_info = elem_info; }
+  virtual void setCurrentElemInfo(const ElemInfo * elem_info);
 
   /**
-   * Set the coordinate system specific volume
+   * Set the coordinate system specific volume, the multiplication with
+   * the transformation factor is done outside of the kernel.
    * @param volume the new coordinate specific volume
    */
   void setCurrentElemVolume(const Real volume) { _current_elem_volume = volume; }
@@ -55,4 +56,7 @@ class LinearFVElementalKernel : public LinearFVKernel
 
   /// The coordinate-specific element volume
   Real _current_elem_volume;
+
+  /// The dof index for the current variable associated with the element
+  dof_id_type _dof_id;
 };
diff --git a/framework/include/linearfvkernels/LinearFVTimeDerivative.h b/framework/include/linearfvkernels/LinearFVTimeDerivative.h
new file mode 100644
index 000000000000..95bc437d9963
--- /dev/null
+++ b/framework/include/linearfvkernels/LinearFVTimeDerivative.h
@@ -0,0 +1,52 @@
+//* 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 "LinearFVElementalKernel.h"
+#include "TimeIntegrator.h"
+
+/**
+ * Kernel that adds contributions from a time derivative term to a linear system
+ * populated using the finite volume method.
+ */
+class LinearFVTimeDerivative : public LinearFVElementalKernel
+{
+public:
+  static InputParameters validParams();
+
+  /**
+   * Class constructor.
+   * @param params The InputParameters for the kernel.
+   */
+  LinearFVTimeDerivative(const InputParameters & params);
+
+  virtual Real computeMatrixContribution() override;
+
+  virtual Real computeRightHandSideContribution() override;
+
+  virtual void setCurrentElemInfo(const ElemInfo * elem_info) override;
+
+protected:
+  /// The functor for the material property multipler
+  const Moose::Functor<Real> & _factor;
+
+  /// The time integrator to use in this kernel, will provide information
+  /// on how many states are required in the history.
+  const TimeIntegrator & _time_integrator;
+
+private:
+  /// Current and older values of the material property multiplier.
+  std::vector<Real> _factor_history;
+
+  /// State args, the args which will help us fetch the different states of
+  /// the material property multiplier. 0th is the current, 1st is old
+  /// 2nd is the older. Might not need all, depends on the time integrator.
+  std::vector<Moose::StateArg> _state_args;
+};
diff --git a/framework/include/timeintegrators/CrankNicolson.h b/framework/include/timeintegrators/CrankNicolson.h
index 90e403707191..1a4eda505670 100644
--- a/framework/include/timeintegrators/CrankNicolson.h
+++ b/framework/include/timeintegrators/CrankNicolson.h
@@ -45,7 +45,7 @@ class CrankNicolson : public TimeIntegrator
 
   virtual Real duDotDuCoeff() const override;
 
-  NumericVector<Number> & _residual_old;
+  NumericVector<Number> * _residual_old;
 };
 
 template <typename T, typename T2>
diff --git a/framework/include/timeintegrators/ExplicitRK2.h b/framework/include/timeintegrators/ExplicitRK2.h
index d4ed2aaed04c..73cdc3eaf5aa 100644
--- a/framework/include/timeintegrators/ExplicitRK2.h
+++ b/framework/include/timeintegrators/ExplicitRK2.h
@@ -79,7 +79,7 @@ class ExplicitRK2 : public TimeIntegrator
   unsigned int _stage;
 
   /// Buffer to store non-time residual from the first stage.
-  NumericVector<Number> & _residual_old;
+  NumericVector<Number> * _residual_old;
 
   /// The older solution
   const NumericVector<Number> & _solution_older;
diff --git a/framework/include/timeintegrators/ExplicitSSPRungeKutta.h b/framework/include/timeintegrators/ExplicitSSPRungeKutta.h
index 5176f9130d90..312c5b292f3c 100644
--- a/framework/include/timeintegrators/ExplicitSSPRungeKutta.h
+++ b/framework/include/timeintegrators/ExplicitSSPRungeKutta.h
@@ -57,9 +57,9 @@ class ExplicitSSPRungeKutta : public ExplicitTimeIntegrator
   std::vector<const NumericVector<Number> *> _solution_stage;
 
   /// Solution vector for intermediate stage
-  NumericVector<Number> & _solution_intermediate_stage;
+  NumericVector<Number> * _solution_intermediate_stage;
   /// Temporary solution vector
-  NumericVector<Number> & _tmp_solution;
+  NumericVector<Number> * _tmp_solution;
   /// Temporary mass-matrix/solution vector product
-  NumericVector<Number> & _tmp_mass_solution_product;
+  NumericVector<Number> * _tmp_mass_solution_product;
 };
diff --git a/framework/include/timeintegrators/ExplicitTVDRK2.h b/framework/include/timeintegrators/ExplicitTVDRK2.h
index b8cd92be8e2a..1fda1fed88d2 100644
--- a/framework/include/timeintegrators/ExplicitTVDRK2.h
+++ b/framework/include/timeintegrators/ExplicitTVDRK2.h
@@ -70,7 +70,7 @@ class ExplicitTVDRK2 : public TimeIntegrator
   unsigned int _stage;
 
   /// Buffer to store non-time residual from the first stage.
-  NumericVector<Number> & _residual_old;
+  NumericVector<Number> * _residual_old;
 
   /// The older solution
   const NumericVector<Number> & _solution_older;
diff --git a/framework/include/timeintegrators/ExplicitTimeIntegrator.h b/framework/include/timeintegrators/ExplicitTimeIntegrator.h
index fd3775bb2cc3..c1eb653c62c9 100644
--- a/framework/include/timeintegrators/ExplicitTimeIntegrator.h
+++ b/framework/include/timeintegrators/ExplicitTimeIntegrator.h
@@ -66,13 +66,13 @@ class ExplicitTimeIntegrator : public TimeIntegrator, public MeshChangedInterfac
   MooseEnum _solve_type;
 
   /// Residual used for the RHS
-  NumericVector<Real> & _explicit_residual;
+  NumericVector<Real> * _explicit_residual;
 
   /// Solution vector for the linear solve
-  NumericVector<Real> & _solution_update;
+  NumericVector<Real> * _solution_update;
 
   /// Diagonal of the lumped mass matrix (and its inversion)
-  NumericVector<Real> & _mass_matrix_diag;
+  NumericVector<Real> * _mass_matrix_diag;
 
   /// Vector of 1's to help with creating the lumped mass matrix
   NumericVector<Real> * _ones;
diff --git a/framework/include/timeintegrators/ImplicitEuler.h b/framework/include/timeintegrators/ImplicitEuler.h
index b71317e5db1b..2e19f8c38e7c 100644
--- a/framework/include/timeintegrators/ImplicitEuler.h
+++ b/framework/include/timeintegrators/ImplicitEuler.h
@@ -30,6 +30,10 @@ class ImplicitEuler : public TimeIntegrator
   virtual void postResidual(NumericVector<Number> & residual) override;
   virtual bool overridesSolve() const override { return false; }
 
+  virtual Real timeDerivativeRHSContribution(const dof_id_type dof_id,
+                                             const std::vector<Real> & factors) const override;
+  virtual Real timeDerivativeMatrixContribution(const Real factor) const override;
+
 protected:
   /**
    * Helper function that actually does the math for computing the time derivative
diff --git a/framework/include/timeintegrators/ImplicitMidpoint.h b/framework/include/timeintegrators/ImplicitMidpoint.h
index bdba146909c3..6869d6afebe4 100644
--- a/framework/include/timeintegrators/ImplicitMidpoint.h
+++ b/framework/include/timeintegrators/ImplicitMidpoint.h
@@ -65,7 +65,7 @@ class ImplicitMidpoint : public TimeIntegrator
   unsigned int _stage;
 
   /// Buffer to store non-time residual from the first stage.
-  NumericVector<Number> & _residual_stage1;
+  NumericVector<Number> * _residual_stage1;
 };
 
 template <typename T, typename T2>
diff --git a/framework/include/timeintegrators/LStableDirk2.h b/framework/include/timeintegrators/LStableDirk2.h
index cc268a50f0bc..e705dbbeb3a9 100644
--- a/framework/include/timeintegrators/LStableDirk2.h
+++ b/framework/include/timeintegrators/LStableDirk2.h
@@ -62,10 +62,10 @@ class LStableDirk2 : public TimeIntegrator
   unsigned int _stage;
 
   //! Buffer to store non-time residual from first stage solve.
-  NumericVector<Number> & _residual_stage1;
+  NumericVector<Number> * _residual_stage1;
 
   //! Buffer to store non-time residual from second stage solve
-  NumericVector<Number> & _residual_stage2;
+  NumericVector<Number> * _residual_stage2;
 
   // The parameter of the method, set at construction time and cannot be changed.
   const Real _alpha;
diff --git a/framework/include/timeintegrators/LinearTimeIntegratorInterface.h b/framework/include/timeintegrators/LinearTimeIntegratorInterface.h
new file mode 100644
index 000000000000..376c65d95b0d
--- /dev/null
+++ b/framework/include/timeintegrators/LinearTimeIntegratorInterface.h
@@ -0,0 +1,52 @@
+//* 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
+
+// Moose includes
+#include "MooseTypes.h"
+
+// Forward declarations
+class SystemBase;
+class LinearSystem;
+
+namespace libMesh
+{
+template <typename T>
+class NumericVector;
+class LinearImplicitSystem;
+} // namespace libMesh
+
+/**
+ * Interface class for routines and member variables for time integrators
+ * relying on linear system assembly method.
+ */
+class LinearTimeIntegratorInterface
+{
+public:
+  LinearTimeIntegratorInterface(SystemBase & system);
+
+  /// The time derivative's contribution to the right hand side of a linear system
+  /// @param dof_id The dof index at which this contribution should be fetched at
+  /// @param factors Multiplicative factor (e.g. a material property) at multiple
+  ///                states (old, older, etc)
+  virtual Real timeDerivativeRHSContribution(dof_id_type dof_id,
+                                             const std::vector<Real> & factors = {}) const;
+
+  /// The time derivative's contribution to the right hand side of a linear system.
+  /// For now, this does not depend of the DoF index, might change in the future.
+  virtual Real timeDerivativeMatrixContribution(const Real factor) const;
+
+protected:
+  /// Pointer to the linear system, can happen that we dont have any
+  LinearSystem * _linear_system;
+
+  /// Nonlinear implicit system, if applicable; otherwise, nullptr
+  libMesh::LinearImplicitSystem * _linear_implicit_system;
+};
diff --git a/framework/include/timeintegrators/NonlinearTimeIntegratorInterface.h b/framework/include/timeintegrators/NonlinearTimeIntegratorInterface.h
new file mode 100644
index 000000000000..9bc56bb7fd1d
--- /dev/null
+++ b/framework/include/timeintegrators/NonlinearTimeIntegratorInterface.h
@@ -0,0 +1,90 @@
+//* 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
+
+// Moose includes
+#include "MooseTypes.h"
+
+// Libmesh includes
+#include "libmesh/enum_parallel_type.h"
+
+// Forward declarations
+class SystemBase;
+class FEProblemBase;
+class NonlinearSystemBase;
+
+namespace libMesh
+{
+template <typename T>
+class NumericVector;
+class NonlinearImplicitSystem;
+} // namespace libMesh
+
+/**
+ * Interface class for routines and member variables for time integrators
+ * relying on Newton's method.
+ */
+class NonlinearTimeIntegratorInterface
+{
+public:
+  NonlinearTimeIntegratorInterface(FEProblemBase & problem, SystemBase & system);
+
+  /**
+   * Callback to the NonLinearTimeIntegratorInterface called immediately after the
+   * residuals are computed in NonlinearSystem::computeResidual().
+   * The residual vector which is passed in to this function should
+   * be filled in by the user with the _Re_time and _Re_non_time
+   * vectors in a way that makes sense for the particular
+   * TimeIntegration method.
+   */
+  virtual void postResidual(NumericVector<Number> & /*residual*/) {}
+
+  /**
+   * Returns the tag for the nodal multiplication factor for the residual calculation of the udot
+   * term.
+   *
+   * By default, this tag will be associated with udot.
+   */
+  TagID uDotFactorTag() const { return _u_dot_factor_tag; }
+  /**
+   * Returns the tag for the nodal multiplication factor for the residual calculation of the udotdot
+   * term.
+   *
+   * By default, this tag will be associated with udotdot.
+   */
+  TagID uDotDotFactorTag() const { return _u_dotdot_factor_tag; }
+
+protected:
+  /// Wrapper around vector addition for nonlinear time integrators. If we don't
+  /// operate on a nonlinear system we don't need to add the vector.
+  /// @param name The name of the vector
+  /// @param project If the vector should be projected
+  /// @param type The parallel distribution of the vetor
+  NumericVector<Number> *
+  addVector(const std::string & name, const bool project, const libMesh::ParallelType type);
+
+  /// Pointer to the nonlinear system, can happen that we dont have any
+  NonlinearSystemBase * _nl;
+
+  /// libMesh nonlinear implicit system, if applicable; otherwise, nullptr
+  libMesh::NonlinearImplicitSystem * _nonlinear_implicit_system;
+
+  /// residual vector for time contributions
+  NumericVector<Number> * _Re_time;
+
+  /// residual vector for non-time contributions
+  NumericVector<Number> * _Re_non_time;
+
+  /// The vector tag for the nodal multiplication factor for the residual calculation of the udot term
+  const TagID _u_dot_factor_tag;
+
+  /// The vector tag for the nodal multiplication factor for the residual calculation of the udotdot term
+  const TagID _u_dotdot_factor_tag;
+};
diff --git a/framework/include/timeintegrators/TimeIntegrator.h b/framework/include/timeintegrators/TimeIntegrator.h
index 3a4c4662005d..e3410a83348a 100644
--- a/framework/include/timeintegrators/TimeIntegrator.h
+++ b/framework/include/timeintegrators/TimeIntegrator.h
@@ -12,6 +12,8 @@
 // MOOSE includes
 #include "MooseObject.h"
 #include "Restartable.h"
+#include "NonlinearTimeIntegratorInterface.h"
+#include "LinearTimeIntegratorInterface.h"
 
 class FEProblemBase;
 class SystemBase;
@@ -21,7 +23,6 @@ namespace libMesh
 {
 template <typename T>
 class NumericVector;
-class NonlinearImplicitSystem;
 } // namespace libMesh
 
 /**
@@ -35,7 +36,10 @@ class NonlinearImplicitSystem;
  * used
  * only by NonlinearSystem (AuxiliarySystem does not produce residual).
  */
-class TimeIntegrator : public MooseObject, public Restartable
+class TimeIntegrator : public MooseObject,
+                       public Restartable,
+                       public NonlinearTimeIntegratorInterface,
+                       public LinearTimeIntegratorInterface
 {
 public:
   static InputParameters validParams();
@@ -69,16 +73,6 @@ class TimeIntegrator : public MooseObject, public Restartable
    */
   virtual void solve();
 
-  /**
-   * Callback to the TimeIntegrator called immediately after the
-   * residuals are computed in NonlinearSystem::computeResidual().
-   * The residual vector which is passed in to this function should
-   * be filled in by the user with the _Re_time and _Re_non_time
-   * vectors in a way that makes sense for the particular
-   * TimeIntegration method.
-   */
-  virtual void postResidual(NumericVector<Number> & /*residual*/) {}
-
   /**
    * Callback to the TimeIntegrator called immediately after
    * TimeIntegrator::solve() (so the name does make sense!).  See
@@ -134,24 +128,15 @@ class TimeIntegrator : public MooseObject, public Restartable
   virtual bool isExplicit() const { return false; }
 
   /**
-   * Returns whether mass matrix is lumped
+   * Return the number of states this requires in a linear
+   * system setting
    */
-  virtual const bool & isLumped() const { return _is_lumped; }
+  virtual unsigned int numStatesRequired() const { return 1; }
 
   /**
-   * Returns the tag for the nodal multiplication factor for the residual calculation of the udot
-   * term.
-   *
-   * By default, this tag will be associated with udot.
-   */
-  TagID uDotFactorTag() const { return _u_dot_factor_tag; }
-  /**
-   * Returns the tag for the nodal multiplication factor for the residual calculation of the udotdot
-   * term.
-   *
-   * By default, this tag will be associated with udotdot.
+   * Returns whether mass matrix is lumped
    */
-  TagID uDotDotFactorTag() const { return _u_dotdot_factor_tag; }
+  virtual const bool & isLumped() const { return _is_lumped; }
 
   /**
    * @returns whether this integrator integrates the given variable
@@ -196,31 +181,32 @@ class TimeIntegrator : public MooseObject, public Restartable
    */
   void computeDuDotDu();
 
+  /// Reference to the problem
   FEProblemBase & _fe_problem;
-  SystemBase & _sys;
-  NonlinearSystemBase & _nl;
-
-  /// Nonlinear implicit system, if applicable; otherwise, nullptr
-  libMesh::NonlinearImplicitSystem * _nonlinear_implicit_system;
 
-  /// residual vector for time contributions
-  NumericVector<Number> & _Re_time;
-  /// residual vector for non-time contributions
-  NumericVector<Number> & _Re_non_time;
+  /// Reference to the system this time integrator operates on
+  SystemBase & _sys;
 
   /// 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
+
+  /// @{
+  /// Solution vectors for different states and variable restrictions
   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;
-  //
+  ///@}
+
+  /// The current time step number
   int & _t_step;
-  //
+
+  /// The current time step size
   Real & _dt;
+
+  /// The previous time step size
   Real & _dt_old;
 
   /// Total number of nonlinear iterations over all stages of the time step
@@ -231,11 +217,6 @@ class TimeIntegrator : public MooseObject, public Restartable
   /// Boolean flag that is set to true if lumped mass matrix is used
   bool _is_lumped;
 
-  /// The vector tag for the nodal multiplication factor for the residual calculation of the udot term
-  const TagID _u_dot_factor_tag;
-  /// The vector tag for the nodal multiplication factor for the residual calculation of the udotdot term
-  const TagID _u_dotdot_factor_tag;
-
   /// Whether the user has requested that the time integrator be applied to a subset of variables
   bool & _var_restriction;
 
diff --git a/framework/src/linearfvkernels/LinearFVElementalKernel.C b/framework/src/linearfvkernels/LinearFVElementalKernel.C
index 8965066af9e3..993ed4a7502c 100644
--- a/framework/src/linearfvkernels/LinearFVElementalKernel.C
+++ b/framework/src/linearfvkernels/LinearFVElementalKernel.C
@@ -29,8 +29,7 @@ LinearFVElementalKernel::addMatrixContribution()
 {
   // These only contribute to the diagonal of the matrix, so we just get
   // the contribution and insert it immediately.
-  const auto dof_id = _current_elem_info->dofIndices()[_sys_num][_var_num];
-  (*_linear_system.matrix).add(dof_id, dof_id, computeMatrixContribution());
+  (*_linear_system.matrix).add(_dof_id, _dof_id, computeMatrixContribution());
 }
 
 void
@@ -38,6 +37,12 @@ LinearFVElementalKernel::addRightHandSideContribution()
 {
   // These only contribute to one entry of the right hand side, so we just get
   // the contribution and insert it immediately.
-  const auto dof_id = _current_elem_info->dofIndices()[_sys_num][_var_num];
-  (*_linear_system.rhs).add(dof_id, computeRightHandSideContribution());
+  (*_linear_system.rhs).add(_dof_id, computeRightHandSideContribution());
+}
+
+void
+LinearFVElementalKernel::setCurrentElemInfo(const ElemInfo * elem_info)
+{
+  _current_elem_info = elem_info;
+  _dof_id = _current_elem_info->dofIndices()[_sys_num][_var_num];
 }
diff --git a/framework/src/linearfvkernels/LinearFVTimeDerivative.C b/framework/src/linearfvkernels/LinearFVTimeDerivative.C
new file mode 100644
index 000000000000..9a95afa2590c
--- /dev/null
+++ b/framework/src/linearfvkernels/LinearFVTimeDerivative.C
@@ -0,0 +1,60 @@
+//* 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 "LinearFVTimeDerivative.h"
+
+registerMooseObject("MooseApp", LinearFVTimeDerivative);
+
+InputParameters
+LinearFVTimeDerivative::validParams()
+{
+  InputParameters params = LinearFVElementalKernel::validParams();
+  params.addClassDescription("Represents the matrix and right hand side contributions of a "
+                             "time derivative term in a partial differential equation.");
+  params.addParam<MooseFunctorName>(
+      "factor", 1.0, "A multiplier on the variable within the time derivative.");
+  return params;
+}
+
+LinearFVTimeDerivative::LinearFVTimeDerivative(const InputParameters & params)
+  : LinearFVElementalKernel(params),
+    _factor(getFunctor<Real>("factor")),
+    _time_integrator(_sys.getTimeIntegrator(_var_num)),
+    _factor_history(_time_integrator.numStatesRequired(), 0.0),
+    _state_args(_time_integrator.numStatesRequired(), determineState())
+{
+  // In case we need older states
+  for (const auto i : index_range(_state_args))
+    _state_args[i] = Moose::StateArg(i, Moose::SolutionIterationType::Time);
+}
+
+Real
+LinearFVTimeDerivative::computeMatrixContribution()
+{
+  const auto elem_arg = makeElemArg(_current_elem_info->elem());
+  return _time_integrator.timeDerivativeMatrixContribution(_factor(elem_arg, _state_args[0])) *
+         _current_elem_volume;
+}
+
+Real
+LinearFVTimeDerivative::computeRightHandSideContribution()
+{
+  return _time_integrator.timeDerivativeRHSContribution(_dof_id, _factor_history) *
+         _current_elem_volume;
+}
+
+void
+LinearFVTimeDerivative::setCurrentElemInfo(const ElemInfo * elem_info)
+{
+  LinearFVElementalKernel::setCurrentElemInfo(elem_info);
+
+  const auto elem_arg = makeElemArg(_current_elem_info->elem());
+  for (const auto i : index_range(_factor_history))
+    _factor_history[i] = _factor(elem_arg, _state_args[i]);
+}
diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C
index 6bb47530c5b5..8e6d293c08db 100644
--- a/framework/src/problems/FEProblemBase.C
+++ b/framework/src/problems/FEProblemBase.C
@@ -6457,13 +6457,13 @@ FEProblemBase::advanceState()
 {
   TIME_SECTION("advanceState", 5, "Advancing State");
 
-  for (auto & nl : _nl)
-    nl->copyOldSolutions();
+  for (auto & sys : _solver_systems)
+    sys->copyOldSolutions();
   _aux->copyOldSolutions();
 
   if (_displaced_problem)
   {
-    for (const auto i : index_range(_nl))
+    for (const auto i : index_range(_solver_systems))
       _displaced_problem->solverSys(i).copyOldSolutions();
     _displaced_problem->auxSys().copyOldSolutions();
   }
diff --git a/framework/src/timeintegrators/AStableDirk4.C b/framework/src/timeintegrators/AStableDirk4.C
index e233631eef00..151f150cbc99 100644
--- a/framework/src/timeintegrators/AStableDirk4.C
+++ b/framework/src/timeintegrators/AStableDirk4.C
@@ -42,7 +42,7 @@ AStableDirk4::AStableDirk4(const InputParameters & parameters)
   {
     std::ostringstream oss;
     oss << "residual_stage" << stage + 1;
-    _stage_residuals[stage] = &(_nl.addVector(oss.str(), false, GHOSTED));
+    _stage_residuals[stage] = addVector(oss.str(), false, GHOSTED);
   }
 
   // Initialize parameters
@@ -148,14 +148,14 @@ AStableDirk4::solve()
       }
 
       // Do the solve
-      _nl.system().solve();
+      _nl->system().solve();
 
       // Update the iteration counts
       _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
       _n_linear_iterations += getNumLinearIterationsLastSolve();
 
       // Abort time step immediately on stage failure - see TimeIntegrator doc page
-      if (!_fe_problem.converged(_nl.number()))
+      if (!_fe_problem.converged(_nl->number()))
         return;
     }
   }
@@ -188,10 +188,10 @@ AStableDirk4::postResidual(NumericVector<Number> & residual)
 
       // Store this stage's non-time residual.  We are calling operator=
       // here, and that calls close().
-      *_stage_residuals[_stage - 1] = _Re_non_time;
+      *_stage_residuals[_stage - 1] = *_Re_non_time;
 
       // Build up the residual for this stage.
-      residual.add(1., _Re_time);
+      residual.add(1., *_Re_time);
       for (unsigned int j = 0; j < _stage; ++j)
         residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
       residual.close();
@@ -206,7 +206,7 @@ AStableDirk4::postResidual(NumericVector<Number> & residual)
       // just do one more stage, but I think handling it separately like
       // this is easier to understand and doesn't create too much code
       // repitition.
-      residual.add(1., _Re_time);
+      residual.add(1., *_Re_time);
       for (unsigned int j = 0; j < 3; ++j)
         residual.add(_b[j], *_stage_residuals[j]);
       residual.close();
diff --git a/framework/src/timeintegrators/ActuallyExplicitEuler.C b/framework/src/timeintegrators/ActuallyExplicitEuler.C
index c6488e3abb61..861b90235ea2 100644
--- a/framework/src/timeintegrators/ActuallyExplicitEuler.C
+++ b/framework/src/timeintegrators/ActuallyExplicitEuler.C
@@ -76,12 +76,12 @@ ActuallyExplicitEuler::solve()
   _nonlinear_implicit_system->update();
 
   // Compute the residual
-  _explicit_residual.zero();
+  _explicit_residual->zero();
   _fe_problem.computeResidual(
-      *_nonlinear_implicit_system->current_local_solution, _explicit_residual, _nl.number());
+      *_nonlinear_implicit_system->current_local_solution, *_explicit_residual, _nl->number());
 
   // Move the residual to the RHS
-  _explicit_residual *= -1.0;
+  *_explicit_residual *= -1.0;
 
   // Compute the mass matrix
   auto & mass_matrix = _nonlinear_implicit_system->get_system_matrix();
@@ -93,14 +93,14 @@ ActuallyExplicitEuler::solve()
   bool converged = performExplicitSolve(mass_matrix);
 
   // Update the solution
-  *_nonlinear_implicit_system->solution = _nl.solutionOld();
-  *_nonlinear_implicit_system->solution += _solution_update;
+  *_nonlinear_implicit_system->solution = _nl->solutionOld();
+  *_nonlinear_implicit_system->solution += *_solution_update;
 
   // Constraints may be solved in an uncoupled way. For example, momentum-balance equations may be
   // solved node-wise and then the solution (e.g. velocities or positions)can be applied to those
   // nodes without solving for such constraints on a system level. This strategy is being used for
   // node-face contact in explicit dynamics.
-  _nl.overwriteNodeFace(*_nonlinear_implicit_system->solution);
+  _nl->overwriteNodeFace(*_nonlinear_implicit_system->solution);
 
   // Enforce contraints on the solution
   DofMap & dof_map = _nonlinear_implicit_system->get_dof_map();
@@ -108,7 +108,7 @@ ActuallyExplicitEuler::solve()
                                       _nonlinear_implicit_system->solution.get());
   _nonlinear_implicit_system->update();
 
-  _nl.setSolution(*_nonlinear_implicit_system->current_local_solution);
+  _nl->setSolution(*_nonlinear_implicit_system->current_local_solution);
 
   _nonlinear_implicit_system->nonlinear_solver->converged = converged;
 }
@@ -116,8 +116,8 @@ ActuallyExplicitEuler::solve()
 void
 ActuallyExplicitEuler::postResidual(NumericVector<Number> & residual)
 {
-  residual += _Re_time;
-  residual += _Re_non_time;
+  residual += *_Re_time;
+  residual += *_Re_non_time;
   residual.close();
 
   // Reset time to the time at which to evaluate nodal BCs, which comes next
diff --git a/framework/src/timeintegrators/BDF2.C b/framework/src/timeintegrators/BDF2.C
index 87f5123df0a1..29487f3558f7 100644
--- a/framework/src/timeintegrators/BDF2.C
+++ b/framework/src/timeintegrators/BDF2.C
@@ -72,8 +72,8 @@ BDF2::computeADTimeDerivatives(ADReal & ad_u_dot,
 void
 BDF2::postResidual(NumericVector<Number> & residual)
 {
-  residual += _Re_time;
-  residual += _Re_non_time;
+  residual += *_Re_time;
+  residual += *_Re_non_time;
   residual.close();
 }
 
diff --git a/framework/src/timeintegrators/CentralDifference.C b/framework/src/timeintegrators/CentralDifference.C
index b6b07393357b..d8d971bf8b05 100644
--- a/framework/src/timeintegrators/CentralDifference.C
+++ b/framework/src/timeintegrators/CentralDifference.C
@@ -61,10 +61,10 @@ CentralDifference::initialSetup()
   ActuallyExplicitEuler::initialSetup();
 
   // _nl here so that we don't create this vector in the aux system time integrator
-  _nl.disassociateVectorFromTag(*_nl.solutionUDot(), _u_dot_factor_tag);
-  _nl.addVector(_u_dot_factor_tag, true, GHOSTED);
-  _nl.disassociateVectorFromTag(*_nl.solutionUDotDot(), _u_dotdot_factor_tag);
-  _nl.addVector(_u_dotdot_factor_tag, true, GHOSTED);
+  _nl->disassociateVectorFromTag(*_nl->solutionUDot(), _u_dot_factor_tag);
+  _nl->addVector(_u_dot_factor_tag, true, GHOSTED);
+  _nl->disassociateVectorFromTag(*_nl->solutionUDotDot(), _u_dotdot_factor_tag);
+  _nl->addVector(_u_dotdot_factor_tag, true, GHOSTED);
 }
 
 void
diff --git a/framework/src/timeintegrators/CrankNicolson.C b/framework/src/timeintegrators/CrankNicolson.C
index 103b4e97e239..64ee40182eb7 100644
--- a/framework/src/timeintegrators/CrankNicolson.C
+++ b/framework/src/timeintegrators/CrankNicolson.C
@@ -22,8 +22,7 @@ CrankNicolson::validParams()
 }
 
 CrankNicolson::CrankNicolson(const InputParameters & parameters)
-  : TimeIntegrator(parameters),
-    _residual_old(_nl.addVector("residual_old", false, libMesh::GHOSTED))
+  : TimeIntegrator(parameters), _residual_old(addVector("residual_old", false, libMesh::GHOSTED))
 {
 }
 
@@ -79,17 +78,20 @@ CrankNicolson::init()
   u_dot.zero();
   computeDuDotDu();
 
-  // compute residual for the initial time step
-  // Note: we can not directly pass _residual_old in computeResidualTag because
-  //       the function will call postResidual, which will cause _residual_old
-  //       to be added on top of itself prohibited by PETSc.
-  //       Objects executed on initial have been executed by FEProblem,
-  //       so we can and should directly call NonlinearSystem residual evaluation.
-  _fe_problem.setCurrentResidualVectorTags({_nl.nonTimeVectorTag()});
-  _nl.computeResidualTag(_nl.RHS(), _nl.nonTimeVectorTag());
-  _fe_problem.clearCurrentResidualVectorTags();
-
-  copyVector(_nl.RHS(), _residual_old);
+  if (_nl)
+  {
+    // compute residual for the initial time step
+    // Note: we can not directly pass _residual_old in computeResidualTag because
+    //       the function will call postResidual, which will cause _residual_old
+    //       to be added on top of itself prohibited by PETSc.
+    //       Objects executed on initial have been executed by FEProblem,
+    //       so we can and should directly call NonlinearSystem residual evaluation.
+    _fe_problem.setCurrentResidualVectorTags({_nl->nonTimeVectorTag()});
+    _nl->computeResidualTag(_nl->RHS(), _nl->nonTimeVectorTag());
+    _fe_problem.clearCurrentResidualVectorTags();
+
+    copyVector(_nl->RHS(), *_residual_old);
+  }
 }
 
 void
@@ -99,7 +101,7 @@ CrankNicolson::postResidual(NumericVector<Number> & residual)
   // and that's probably a good idea with earlier versions too, but
   // we don't always get here with _Re_time closed.
   std::vector<unsigned char> inputs_closed = {
-      _Re_time.closed(), _Re_non_time.closed(), _residual_old.closed()};
+      _Re_time->closed(), _Re_non_time->closed(), _residual_old->closed()};
 
   // We might have done work on one processor but not all processors,
   // so we have to sync our closed() checks.  Congrats to the BISON
@@ -107,38 +109,38 @@ CrankNicolson::postResidual(NumericVector<Number> & residual)
   comm().min(inputs_closed);
 
   if (!inputs_closed[0])
-    _Re_time.close();
+    _Re_time->close();
   if (!inputs_closed[1])
-    _Re_non_time.close();
+    _Re_non_time->close();
   if (!inputs_closed[2])
-    _residual_old.close();
+    _residual_old->close();
 
   if (!_var_restriction)
   {
-    residual += _Re_time;
-    residual += _Re_non_time;
-    residual += _residual_old;
+    residual += *_Re_time;
+    residual += *_Re_non_time;
+    residual += *_residual_old;
   }
   else
   {
     auto residual_sub = residual.get_subvector(_local_indices);
-    auto re_time_sub = _Re_time.get_subvector(_local_indices);
-    auto re_non_time_sub = _Re_non_time.get_subvector(_local_indices);
-    auto residual_old_sub = _residual_old.get_subvector(_local_indices);
+    auto re_time_sub = _Re_time->get_subvector(_local_indices);
+    auto re_non_time_sub = _Re_non_time->get_subvector(_local_indices);
+    auto residual_old_sub = _residual_old->get_subvector(_local_indices);
     *residual_sub += *re_time_sub;
     *residual_sub += *re_non_time_sub;
     *residual_sub += *residual_old_sub;
     residual.restore_subvector(std::move(residual_sub), _local_indices);
-    _Re_time.restore_subvector(std::move(re_time_sub), _local_indices);
-    _Re_non_time.restore_subvector(std::move(re_non_time_sub), _local_indices);
-    _residual_old.restore_subvector(std::move(residual_old_sub), _local_indices);
+    _Re_time->restore_subvector(std::move(re_time_sub), _local_indices);
+    _Re_non_time->restore_subvector(std::move(re_non_time_sub), _local_indices);
+    _residual_old->restore_subvector(std::move(residual_old_sub), _local_indices);
   }
 }
 
 void
 CrankNicolson::postStep()
 {
-  copyVector(_Re_non_time, _residual_old);
+  copyVector(*_Re_non_time, *_residual_old);
 }
 
 Real
diff --git a/framework/src/timeintegrators/ExplicitEuler.C b/framework/src/timeintegrators/ExplicitEuler.C
index 26f0cd76d98e..76ac3d149d59 100644
--- a/framework/src/timeintegrators/ExplicitEuler.C
+++ b/framework/src/timeintegrators/ExplicitEuler.C
@@ -57,7 +57,7 @@ ExplicitEuler::computeADTimeDerivatives(ADReal & ad_u_dot,
 void
 ExplicitEuler::postResidual(NumericVector<Number> & residual)
 {
-  residual += _Re_time;
-  residual += _Re_non_time;
+  residual += *_Re_time;
+  residual += *_Re_non_time;
   residual.close();
 }
diff --git a/framework/src/timeintegrators/ExplicitRK2.C b/framework/src/timeintegrators/ExplicitRK2.C
index bb5ed9e1920e..76107e1f4618 100644
--- a/framework/src/timeintegrators/ExplicitRK2.C
+++ b/framework/src/timeintegrators/ExplicitRK2.C
@@ -25,7 +25,7 @@ ExplicitRK2::validParams()
 ExplicitRK2::ExplicitRK2(const InputParameters & parameters)
   : TimeIntegrator(parameters),
     _stage(1),
-    _residual_old(_nl.addVector("residual_old", false, GHOSTED)),
+    _residual_old(addVector("residual_old", false, GHOSTED)),
     _solution_older(_sys.solutionState(2))
 {
   mooseInfo("ExplicitRK2-derived TimeIntegrators (ExplicitMidpoint, Heun, Ralston) and other "
@@ -86,12 +86,12 @@ ExplicitRK2::solve()
   _stage = 2;
   _fe_problem.timeOld() = time_old;
   _fe_problem.time() = time_stage2;
-  _nl.system().solve();
+  _nl->system().solve();
   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
   _n_linear_iterations += getNumLinearIterationsLastSolve();
 
   // Abort time step immediately on stage failure - see TimeIntegrator doc page
-  if (!_fe_problem.converged(_nl.number()))
+  if (!_fe_problem.converged(_nl->number()))
     return;
 
   // Advance solutions old->older, current->old.  Also moves Material
@@ -105,7 +105,7 @@ ExplicitRK2::solve()
   _stage = 3;
   _fe_problem.timeOld() = time_stage2;
   _fe_problem.time() = time_new;
-  _nl.system().solve();
+  _nl->system().solve();
   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
   _n_linear_iterations += getNumLinearIterationsLastSolve();
 
@@ -136,11 +136,11 @@ ExplicitRK2::postResidual(NumericVector<Number> & residual)
     // .) The minus signs are "baked in" to the non-time residuals, so
     //    they do not appear here.
     // .) The current non-time residual is saved for the next stage.
-    _residual_old = _Re_non_time;
-    _residual_old.close();
+    *_residual_old = *_Re_non_time;
+    _residual_old->close();
 
-    residual.add(1., _Re_time);
-    residual.add(a(), _residual_old);
+    residual.add(1., *_Re_time);
+    residual.add(a(), *_residual_old);
     residual.close();
   }
   else if (_stage == 3)
@@ -156,9 +156,9 @@ ExplicitRK2::postResidual(NumericVector<Number> & residual)
     //    residuals, so it does not appear here.
     // .) Although this is an update step, we have to do a "solve"
     //    using the mass matrix.
-    residual.add(1., _Re_time);
-    residual.add(b1(), _residual_old);
-    residual.add(b2(), _Re_non_time);
+    residual.add(1., *_Re_time);
+    residual.add(b1(), *_residual_old);
+    residual.add(b2(), *_Re_non_time);
     residual.close();
   }
   else
diff --git a/framework/src/timeintegrators/ExplicitSSPRungeKutta.C b/framework/src/timeintegrators/ExplicitSSPRungeKutta.C
index c1da27e2c998..3c152f5058ed 100644
--- a/framework/src/timeintegrators/ExplicitSSPRungeKutta.C
+++ b/framework/src/timeintegrators/ExplicitSSPRungeKutta.C
@@ -32,13 +32,11 @@ ExplicitSSPRungeKutta::validParams()
 
 ExplicitSSPRungeKutta::ExplicitSSPRungeKutta(const InputParameters & parameters)
   : ExplicitTimeIntegrator(parameters),
-
     _order(getParam<MooseEnum>("order")),
-
     _stage(0),
-    _solution_intermediate_stage(_nl.addVector("solution_intermediate_stage", false, GHOSTED)),
-    _tmp_solution(_nl.addVector("tmp_solution", false, GHOSTED)),
-    _tmp_mass_solution_product(_nl.addVector("tmp_mass_solution_product", false, GHOSTED))
+    _solution_intermediate_stage(addVector("solution_intermediate_stage", false, GHOSTED)),
+    _tmp_solution(addVector("tmp_solution", false, GHOSTED)),
+    _tmp_mass_solution_product(addVector("tmp_mass_solution_product", false, GHOSTED))
 {
   // For SSPRK methods up to order 3, the number of stages equals the order
   _n_stages = _order;
@@ -125,9 +123,9 @@ ExplicitSSPRungeKutta::solve()
     else
     {
       // Else must be the intermediate stage of the 3-stage method
-      _solution_intermediate_stage = *_solution;
-      _solution_intermediate_stage.close();
-      _solution_stage[_stage] = &_solution_intermediate_stage;
+      *_solution_intermediate_stage = *_solution;
+      _solution_intermediate_stage->close();
+      _solution_stage[_stage] = _solution_intermediate_stage;
     }
 
     // Set stage time for residual evaluation
@@ -156,12 +154,12 @@ ExplicitSSPRungeKutta::solveStage()
       *_nonlinear_implicit_system->current_local_solution, mass_matrix, _Ke_time_tag);
 
   // Compute RHS vector using previous stage solution in steady-state residual
-  _explicit_residual.zero();
+  _explicit_residual->zero();
   _fe_problem.computeResidual(
-      *_nonlinear_implicit_system->current_local_solution, _explicit_residual, _nl.number());
+      *_nonlinear_implicit_system->current_local_solution, *_explicit_residual, _nl->number());
 
   // Move the residual to the RHS
-  _explicit_residual *= -1.0;
+  *_explicit_residual *= -1.0;
 
   // Perform the linear solve
   bool converged = performExplicitSolve(mass_matrix);
@@ -169,7 +167,7 @@ ExplicitSSPRungeKutta::solveStage()
   // Update the solution: u^(s) = u^(s-1) + du^(s)
   (*_nonlinear_implicit_system->solution).zero();
   *_nonlinear_implicit_system->solution = *(_solution_stage[_stage]);
-  *_nonlinear_implicit_system->solution += _solution_update;
+  *_nonlinear_implicit_system->solution += *_solution_update;
 
   // Enforce contraints on the solution
   DofMap & dof_map = _nonlinear_implicit_system->get_dof_map();
@@ -177,7 +175,7 @@ ExplicitSSPRungeKutta::solveStage()
                                       _nonlinear_implicit_system->solution.get());
   _nonlinear_implicit_system->update();
 
-  _nl.setSolution(*_nonlinear_implicit_system->current_local_solution);
+  _nl->setSolution(*_nonlinear_implicit_system->current_local_solution);
 
   _nonlinear_implicit_system->nonlinear_solver->converged = converged;
 
@@ -188,21 +186,21 @@ void
 ExplicitSSPRungeKutta::postResidual(NumericVector<Number> & residual)
 {
   // The time residual is not included in the steady-state residual
-  residual += _Re_non_time;
+  residual += *_Re_non_time;
 
   // Compute \sum_{k=0}^{s-1} a_{s,k} u^(k) - u^(s-1)
-  _tmp_solution.zero();
+  _tmp_solution->zero();
   for (unsigned int k = 0; k <= _stage; k++)
-    _tmp_solution.add(_a[_stage][k], *(_solution_stage[k]));
-  _tmp_solution.add(-1.0, *(_solution_stage[_stage]));
-  _tmp_solution.close();
+    _tmp_solution->add(_a[_stage][k], *(_solution_stage[k]));
+  _tmp_solution->add(-1.0, *(_solution_stage[_stage]));
+  _tmp_solution->close();
 
   // Perform mass matrix product with the above vector
   auto & mass_matrix = _nonlinear_implicit_system->get_system_matrix();
-  mass_matrix.vector_mult(_tmp_mass_solution_product, _tmp_solution);
+  mass_matrix.vector_mult(*_tmp_mass_solution_product, *_tmp_solution);
 
   // Finish computing residual vector (before modification by nodal BCs)
-  residual -= _tmp_mass_solution_product;
+  residual -= *_tmp_mass_solution_product;
 
   residual.close();
 
diff --git a/framework/src/timeintegrators/ExplicitTVDRK2.C b/framework/src/timeintegrators/ExplicitTVDRK2.C
index 61d4da4a66bb..336486955685 100644
--- a/framework/src/timeintegrators/ExplicitTVDRK2.C
+++ b/framework/src/timeintegrators/ExplicitTVDRK2.C
@@ -26,7 +26,7 @@ ExplicitTVDRK2::validParams()
 ExplicitTVDRK2::ExplicitTVDRK2(const InputParameters & parameters)
   : TimeIntegrator(parameters),
     _stage(1),
-    _residual_old(_nl.addVector("residual_old", false, libMesh::GHOSTED)),
+    _residual_old(addVector("residual_old", false, libMesh::GHOSTED)),
     _solution_older(_sys.solutionState(2))
 {
   mooseInfo("ExplicitTVDRK2 and other multistage TimeIntegrators are known not to work with "
@@ -88,12 +88,12 @@ ExplicitTVDRK2::solve()
   _stage = 2;
   _fe_problem.timeOld() = time_old;
   _fe_problem.time() = time_stage2;
-  _nl.system().solve();
+  _nl->system().solve();
   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
   _n_linear_iterations += getNumLinearIterationsLastSolve();
 
   // Abort time step immediately on stage failure - see TimeIntegrator doc page
-  if (!_fe_problem.converged(_nl.number()))
+  if (!_fe_problem.converged(_nl->number()))
     return;
 
   // Advance solutions old->older, current->old.  Also moves Material
@@ -107,7 +107,7 @@ ExplicitTVDRK2::solve()
   _stage = 3;
   _fe_problem.timeOld() = time_stage2;
   _fe_problem.time() = time_new;
-  _nl.system().solve();
+  _nl->system().solve();
   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
   _n_linear_iterations += getNumLinearIterationsLastSolve();
 
@@ -138,11 +138,11 @@ ExplicitTVDRK2::postResidual(NumericVector<Number> & residual)
     // .) The minus signs are "baked in" to the non-time residuals, so
     //    they do not appear here.
     // .) The current non-time residual is saved for the next stage.
-    _residual_old = _Re_non_time;
-    _residual_old.close();
+    *_residual_old = *_Re_non_time;
+    _residual_old->close();
 
-    residual.add(1.0, _Re_time);
-    residual.add(1.0, _residual_old);
+    residual.add(1.0, *_Re_time);
+    residual.add(1.0, *_residual_old);
     residual.close();
   }
   else if (_stage == 3)
@@ -158,8 +158,8 @@ ExplicitTVDRK2::postResidual(NumericVector<Number> & residual)
     //    residuals, so it does not appear here.
     // .) Although this is an update step, we have to do a "solve"
     //    using the mass matrix.
-    residual.add(1.0, _Re_time);
-    residual.add(0.5, _Re_non_time);
+    residual.add(1.0, *_Re_time);
+    residual.add(0.5, *_Re_non_time);
     residual.close();
   }
   else
diff --git a/framework/src/timeintegrators/ExplicitTimeIntegrator.C b/framework/src/timeintegrators/ExplicitTimeIntegrator.C
index 53a2a9844fa7..616ddffbfcc2 100644
--- a/framework/src/timeintegrators/ExplicitTimeIntegrator.C
+++ b/framework/src/timeintegrators/ExplicitTimeIntegrator.C
@@ -41,9 +41,9 @@ ExplicitTimeIntegrator::ExplicitTimeIntegrator(const InputParameters & parameter
     MeshChangedInterface(parameters),
 
     _solve_type(getParam<MooseEnum>("solve_type")),
-    _explicit_residual(_nl.addVector("explicit_residual", false, PARALLEL)),
-    _solution_update(_nl.addVector("solution_update", true, PARALLEL)),
-    _mass_matrix_diag(_nl.addVector("mass_matrix_diag", false, PARALLEL))
+    _explicit_residual(addVector("explicit_residual", false, PARALLEL)),
+    _solution_update(addVector("solution_update", true, PARALLEL)),
+    _mass_matrix_diag(addVector("mass_matrix_diag", false, PARALLEL))
 {
   _Ke_time_tag = _fe_problem.getMatrixTagID("TIME");
 
@@ -52,7 +52,7 @@ ExplicitTimeIntegrator::ExplicitTimeIntegrator(const InputParameters & parameter
   _fe_problem.solverParams()._type = Moose::ST_LINEAR;
 
   if (_solve_type == LUMPED || _solve_type == LUMP_PRECONDITIONED)
-    _ones = &_nl.addVector("ones", false, PARALLEL);
+    _ones = addVector("ones", false, PARALLEL);
 }
 
 void
@@ -86,7 +86,7 @@ ExplicitTimeIntegrator::meshChanged()
 
   if (_solve_type == LUMP_PRECONDITIONED)
   {
-    _preconditioner = std::make_unique<LumpedPreconditioner>(_mass_matrix_diag);
+    _preconditioner = std::make_unique<LumpedPreconditioner>(*_mass_matrix_diag);
     _linear_solver->attach_preconditioner(_preconditioner.get());
     _linear_solver->init();
   }
@@ -113,16 +113,16 @@ ExplicitTimeIntegrator::performExplicitSolve(SparseMatrix<Number> & mass_matrix)
       // Computes the sum of each row (lumping)
       // Note: This is actually how PETSc does it
       // It's not "perfectly optimal" - but it will be fast (and universal)
-      mass_matrix.vector_mult(_mass_matrix_diag, *_ones);
+      mass_matrix.vector_mult(*_mass_matrix_diag, *_ones);
 
       // "Invert" the diagonal mass matrix
-      _mass_matrix_diag.reciprocal();
+      _mass_matrix_diag->reciprocal();
 
       // Multiply the inversion by the RHS
-      _solution_update.pointwise_mult(_mass_matrix_diag, _explicit_residual);
+      _solution_update->pointwise_mult(*_mass_matrix_diag, *_explicit_residual);
 
       // Check for convergence by seeing if there is a nan or inf
-      auto sum = _solution_update.sum();
+      auto sum = _solution_update->sum();
       converged = std::isfinite(sum);
 
       // The linear iteration count remains zero
@@ -132,8 +132,8 @@ ExplicitTimeIntegrator::performExplicitSolve(SparseMatrix<Number> & mass_matrix)
     }
     case LUMP_PRECONDITIONED:
     {
-      mass_matrix.vector_mult(_mass_matrix_diag, *_ones);
-      _mass_matrix_diag.reciprocal();
+      mass_matrix.vector_mult(*_mass_matrix_diag, *_ones);
+      _mass_matrix_diag->reciprocal();
 
       converged = solveLinearSystem(mass_matrix);
 
@@ -153,8 +153,8 @@ ExplicitTimeIntegrator::solveLinearSystem(SparseMatrix<Number> & mass_matrix)
 
   const auto num_its_and_final_tol =
       _linear_solver->solve(mass_matrix,
-                            _solution_update,
-                            _explicit_residual,
+                            *_solution_update,
+                            *_explicit_residual,
                             es.parameters.get<Real>("linear solver tolerance"),
                             es.parameters.get<unsigned int>("linear solver maximum iterations"));
 
diff --git a/framework/src/timeintegrators/ImplicitEuler.C b/framework/src/timeintegrators/ImplicitEuler.C
index 5a2a6a447781..a9f22b9e3d95 100644
--- a/framework/src/timeintegrators/ImplicitEuler.C
+++ b/framework/src/timeintegrators/ImplicitEuler.C
@@ -65,19 +65,34 @@ ImplicitEuler::postResidual(NumericVector<Number> & residual)
 {
   if (!_var_restriction)
   {
-    residual += _Re_time;
-    residual += _Re_non_time;
+    residual += *_Re_time;
+    residual += *_Re_non_time;
     residual.close();
   }
   else
   {
     auto residual_sub = residual.get_subvector(_local_indices);
-    auto re_time_sub = _Re_time.get_subvector(_local_indices);
-    auto re_non_time_sub = _Re_non_time.get_subvector(_local_indices);
+    auto re_time_sub = _Re_time->get_subvector(_local_indices);
+    auto re_non_time_sub = _Re_non_time->get_subvector(_local_indices);
     *residual_sub += *re_time_sub;
     *residual_sub += *re_non_time_sub;
     residual.restore_subvector(std::move(residual_sub), _local_indices);
-    _Re_time.restore_subvector(std::move(re_time_sub), _local_indices);
-    _Re_non_time.restore_subvector(std::move(re_non_time_sub), _local_indices);
+    _Re_time->restore_subvector(std::move(re_time_sub), _local_indices);
+    _Re_non_time->restore_subvector(std::move(re_non_time_sub), _local_indices);
   }
 }
+
+Real
+ImplicitEuler::timeDerivativeRHSContribution(dof_id_type dof_id,
+                                             const std::vector<Real> & factors) const
+{
+  mooseAssert(factors.size() == numStatesRequired(),
+              "Either too many or too few states are given!");
+  return factors[0] * _solution_old(dof_id) / _dt;
+}
+
+Real
+ImplicitEuler::timeDerivativeMatrixContribution(const Real factor) const
+{
+  return factor / _dt;
+}
diff --git a/framework/src/timeintegrators/ImplicitMidpoint.C b/framework/src/timeintegrators/ImplicitMidpoint.C
index 5bde2899fdec..e745b5ce0e03 100644
--- a/framework/src/timeintegrators/ImplicitMidpoint.C
+++ b/framework/src/timeintegrators/ImplicitMidpoint.C
@@ -25,7 +25,7 @@ ImplicitMidpoint::validParams()
 ImplicitMidpoint::ImplicitMidpoint(const InputParameters & parameters)
   : TimeIntegrator(parameters),
     _stage(1),
-    _residual_stage1(_nl.addVector("residual_stage1", false, libMesh::GHOSTED))
+    _residual_stage1(addVector("residual_stage1", false, libMesh::GHOSTED))
 {
   mooseInfo("ImplicitMidpoint and other multistage TimeIntegrators are known not to work with "
             "Materials/AuxKernels that accumulate 'state' and should be used with caution.");
@@ -72,12 +72,12 @@ ImplicitMidpoint::solve()
   _console << "1st stage" << std::endl;
   _stage = 1;
   _fe_problem.time() = time_half;
-  _nl.system().solve();
+  _nl->system().solve();
   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
   _n_linear_iterations += getNumLinearIterationsLastSolve();
 
   // Abort time step immediately on stage failure - see TimeIntegrator doc page
-  if (!_fe_problem.converged(_nl.number()))
+  if (!_fe_problem.converged(_nl->number()))
     return;
 
   // Compute second stage
@@ -85,7 +85,7 @@ ImplicitMidpoint::solve()
   _console << "2nd stage" << std::endl;
   _stage = 2;
   _fe_problem.time() = time_new;
-  _nl.system().solve();
+  _nl->system().solve();
   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
   _n_linear_iterations += getNumLinearIterationsLastSolve();
 }
@@ -103,11 +103,11 @@ ImplicitMidpoint::postResidual(NumericVector<Number> & residual)
     // .) M is the mass matrix
     // .) f(t_n + dt/2, Y_1) is saved in _residual_stage1
     // .) The minus sign is baked in to the non-time residuals, so it does not appear here.
-    _residual_stage1 = _Re_non_time;
-    _residual_stage1.close();
+    *_residual_stage1 = *_Re_non_time;
+    _residual_stage1->close();
 
-    residual.add(1., _Re_time);
-    residual.add(0.5, _Re_non_time);
+    residual.add(1., *_Re_time);
+    residual.add(0.5, *_Re_non_time);
     residual.close();
   }
   else if (_stage == 2)
@@ -123,8 +123,8 @@ ImplicitMidpoint::postResidual(NumericVector<Number> & residual)
     //    been saved as _residual_stage1.
     // .) The minus signs are "baked in" to the non-time residuals, so
     //    they do not appear here.
-    residual.add(1., _Re_time);
-    residual.add(1., _residual_stage1);
+    residual.add(1., *_Re_time);
+    residual.add(1., *_residual_stage1);
     residual.close();
   }
   else
diff --git a/framework/src/timeintegrators/LStableDirk2.C b/framework/src/timeintegrators/LStableDirk2.C
index f84f1e4aae89..536f083c6fff 100644
--- a/framework/src/timeintegrators/LStableDirk2.C
+++ b/framework/src/timeintegrators/LStableDirk2.C
@@ -28,8 +28,8 @@ LStableDirk2::validParams()
 LStableDirk2::LStableDirk2(const InputParameters & parameters)
   : TimeIntegrator(parameters),
     _stage(1),
-    _residual_stage1(_nl.addVector("residual_stage1", false, GHOSTED)),
-    _residual_stage2(_nl.addVector("residual_stage2", false, GHOSTED)),
+    _residual_stage1(addVector("residual_stage1", false, GHOSTED)),
+    _residual_stage2(addVector("residual_stage2", false, GHOSTED)),
     _alpha(1. - 0.5 * std::sqrt(2))
 {
   mooseInfo("LStableDirk2 and other multistage TimeIntegrators are known not to work with "
@@ -82,13 +82,13 @@ LStableDirk2::solve()
   _console << "1st stage" << std::endl;
   _stage = 1;
   _fe_problem.time() = time_stage1;
-  _nl.system().solve();
-  _nl.destroyColoring();
+  _nl->system().solve();
+  _nl->destroyColoring();
   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
   _n_linear_iterations += getNumLinearIterationsLastSolve();
 
   // Abort time step immediately on stage failure - see TimeIntegrator doc page
-  if (!_fe_problem.converged(_nl.number()))
+  if (!_fe_problem.converged(_nl->number()))
     return;
 
   // Compute second stage
@@ -97,8 +97,8 @@ LStableDirk2::solve()
   _stage = 2;
   _fe_problem.timeOld() = time_stage1;
   _fe_problem.time() = time_new;
-  _nl.potentiallySetupFiniteDifferencing();
-  _nl.system().solve();
+  _nl->potentiallySetupFiniteDifferencing();
+  _nl->system().solve();
   _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
   _n_linear_iterations += getNumLinearIterationsLastSolve();
 
@@ -122,11 +122,11 @@ LStableDirk2::postResidual(NumericVector<Number> & residual)
     // .) (Y_1 - y_n)/dt corresponds to the residual of the time kernels.
     // .) The minus sign in front of alpha is already "baked in" to
     //    the non-time residuals, so it does not appear here.
-    _residual_stage1 = _Re_non_time;
-    _residual_stage1.close();
+    *_residual_stage1 = *_Re_non_time;
+    _residual_stage1->close();
 
-    residual.add(1., _Re_time);
-    residual.add(_alpha, _residual_stage1);
+    residual.add(1., *_Re_time);
+    residual.add(_alpha, *_residual_stage1);
     residual.close();
   }
   else if (_stage == 2)
@@ -143,12 +143,12 @@ LStableDirk2::postResidual(NumericVector<Number> & residual)
     //    residuals, so they do not appear here.
     //
     // The solution at the end of stage 2, i.e. Y_2, is also the final solution.
-    _residual_stage2 = _Re_non_time;
-    _residual_stage2.close();
+    *_residual_stage2 = *_Re_non_time;
+    _residual_stage2->close();
 
-    residual.add(1., _Re_time);
-    residual.add(1. - _alpha, _residual_stage1);
-    residual.add(_alpha, _residual_stage2);
+    residual.add(1., *_Re_time);
+    residual.add(1. - _alpha, *_residual_stage1);
+    residual.add(_alpha, *_residual_stage2);
     residual.close();
   }
   else
diff --git a/framework/src/timeintegrators/LStableDirk3.C b/framework/src/timeintegrators/LStableDirk3.C
index 4b6771778c6a..4449399e7bdd 100644
--- a/framework/src/timeintegrators/LStableDirk3.C
+++ b/framework/src/timeintegrators/LStableDirk3.C
@@ -38,7 +38,7 @@ LStableDirk3::LStableDirk3(const InputParameters & parameters)
   {
     std::ostringstream oss;
     oss << "residual_stage" << stage + 1;
-    _stage_residuals[stage] = &(_nl.addVector(oss.str(), false, GHOSTED));
+    _stage_residuals[stage] = addVector(oss.str(), false, GHOSTED);
   }
 
   // Initialize parameters
@@ -108,20 +108,20 @@ LStableDirk3::solve()
 
     // If we previously used coloring, destroy the old object so it doesn't leak when we allocate a
     // new object in the following lines
-    _nl.destroyColoring();
+    _nl->destroyColoring();
 
     // Potentially setup finite differencing contexts for the solve
-    _nl.potentiallySetupFiniteDifferencing();
+    _nl->potentiallySetupFiniteDifferencing();
 
     // Do the solve
-    _nl.system().solve();
+    _nl->system().solve();
 
     // Update the iteration counts
     _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
     _n_linear_iterations += getNumLinearIterationsLastSolve();
 
     // Abort time step immediately on stage failure - see TimeIntegrator doc page
-    if (!_fe_problem.converged(_nl.number()))
+    if (!_fe_problem.converged(_nl->number()))
       return;
   }
 }
@@ -147,10 +147,10 @@ LStableDirk3::postResidual(NumericVector<Number> & residual)
 
   // Store this stage's non-time residual.  We are calling operator=
   // here, and that calls close().
-  *_stage_residuals[_stage - 1] = _Re_non_time;
+  *_stage_residuals[_stage - 1] = *_Re_non_time;
 
   // Build up the residual for this stage.
-  residual.add(1., _Re_time);
+  residual.add(1., *_Re_time);
   for (unsigned int j = 0; j < _stage; ++j)
     residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
   residual.close();
diff --git a/framework/src/timeintegrators/LStableDirk4.C b/framework/src/timeintegrators/LStableDirk4.C
index 202a993f6ae3..159653076d74 100644
--- a/framework/src/timeintegrators/LStableDirk4.C
+++ b/framework/src/timeintegrators/LStableDirk4.C
@@ -46,7 +46,7 @@ LStableDirk4::LStableDirk4(const InputParameters & parameters)
   {
     std::ostringstream oss;
     oss << "residual_stage" << stage + 1;
-    _stage_residuals[stage] = &(_nl.addVector(oss.str(), false, GHOSTED));
+    _stage_residuals[stage] = addVector(oss.str(), false, GHOSTED);
   }
 }
 
@@ -104,20 +104,20 @@ LStableDirk4::solve()
 
     // If we previously used coloring, destroy the old object so it doesn't leak when we allocate a
     // new object in the following lines
-    _nl.destroyColoring();
+    _nl->destroyColoring();
 
     // Potentially setup finite differencing contexts for the solve
-    _nl.potentiallySetupFiniteDifferencing();
+    _nl->potentiallySetupFiniteDifferencing();
 
     // Do the solve
-    _nl.system().solve();
+    _nl->system().solve();
 
     // Update the iteration counts
     _n_nonlinear_iterations += getNumNonlinearIterationsLastSolve();
     _n_linear_iterations += getNumLinearIterationsLastSolve();
 
     // Abort time step immediately on stage failure - see TimeIntegrator doc page
-    if (!_fe_problem.converged(_nl.number()))
+    if (!_fe_problem.converged(_nl->number()))
       return;
   }
 }
@@ -146,10 +146,10 @@ LStableDirk4::postResidual(NumericVector<Number> & residual)
 
   // Store this stage's non-time residual.  We are calling operator=
   // here, and that calls close().
-  *_stage_residuals[_stage - 1] = _Re_non_time;
+  *_stage_residuals[_stage - 1] = *_Re_non_time;
 
   // Build up the residual for this stage.
-  residual.add(1., _Re_time);
+  residual.add(1., *_Re_time);
   for (unsigned int j = 0; j < _stage; ++j)
     residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
   residual.close();
diff --git a/framework/src/timeintegrators/LinearTimeIntegratorInterface.C b/framework/src/timeintegrators/LinearTimeIntegratorInterface.C
new file mode 100644
index 000000000000..da94a504837e
--- /dev/null
+++ b/framework/src/timeintegrators/LinearTimeIntegratorInterface.C
@@ -0,0 +1,38 @@
+//* 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 "LinearTimeIntegratorInterface.h"
+#include "FEProblem.h"
+#include "LinearSystem.h"
+
+#include "libmesh/linear_implicit_system.h"
+
+LinearTimeIntegratorInterface::LinearTimeIntegratorInterface(SystemBase & system)
+  : _linear_system(dynamic_cast<LinearSystem *>(&system)),
+    _linear_implicit_system(
+        _linear_system ? dynamic_cast<LinearImplicitSystem *>(&_linear_system->system()) : nullptr)
+{
+}
+
+Real
+LinearTimeIntegratorInterface::timeDerivativeRHSContribution(
+    const dof_id_type /*dof_id*/, const std::vector<Real> & /*factors*/) const
+{
+  mooseError("The time derivative right hand side contribution has not been implemented yet",
+             _linear_system ? " for time integrator of system " + _linear_system->name() : "",
+             "!");
+}
+
+Real
+LinearTimeIntegratorInterface::timeDerivativeMatrixContribution(const Real /*factor*/) const
+{
+  mooseError("The time derivative matrix contribution has not been implemented yet",
+             _linear_system ? " for time integrator of system " + _linear_system->name() : "",
+             "!");
+}
diff --git a/framework/src/timeintegrators/NewmarkBeta.C b/framework/src/timeintegrators/NewmarkBeta.C
index 4ff7475fe254..d02b8302619f 100644
--- a/framework/src/timeintegrators/NewmarkBeta.C
+++ b/framework/src/timeintegrators/NewmarkBeta.C
@@ -114,8 +114,8 @@ NewmarkBeta::computeADTimeDerivatives(ADReal & ad_u_dot,
 void
 NewmarkBeta::postResidual(NumericVector<Number> & residual)
 {
-  residual += _Re_time;
-  residual += _Re_non_time;
+  residual += *_Re_time;
+  residual += *_Re_non_time;
   residual.close();
 }
 
diff --git a/framework/src/timeintegrators/NonlinearTimeIntegratorInterface.C b/framework/src/timeintegrators/NonlinearTimeIntegratorInterface.C
new file mode 100644
index 000000000000..7784e6c98df2
--- /dev/null
+++ b/framework/src/timeintegrators/NonlinearTimeIntegratorInterface.C
@@ -0,0 +1,39 @@
+//* 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 "NonlinearTimeIntegratorInterface.h"
+#include "FEProblem.h"
+#include "NonlinearSystemBase.h"
+
+#include "libmesh/nonlinear_implicit_system.h"
+#include "libmesh/nonlinear_solver.h"
+#include "libmesh/dof_map.h"
+
+NonlinearTimeIntegratorInterface::NonlinearTimeIntegratorInterface(FEProblemBase & problem,
+                                                                   SystemBase & system)
+  : _nl(dynamic_cast<NonlinearSystemBase *>(&system)),
+    _nonlinear_implicit_system(_nl ? dynamic_cast<NonlinearImplicitSystem *>(&_nl->system())
+                                   : nullptr),
+    _Re_time(_nl ? &_nl->getResidualTimeVector() : nullptr),
+    _Re_non_time(_nl ? &_nl->getResidualNonTimeVector() : nullptr),
+    _u_dot_factor_tag(problem.addVectorTag("u_dot_factor", Moose::VECTOR_TAG_SOLUTION)),
+    _u_dotdot_factor_tag(problem.addVectorTag("u_dotdot_factor", Moose::VECTOR_TAG_SOLUTION))
+{
+}
+
+NumericVector<Number> *
+NonlinearTimeIntegratorInterface::addVector(const std::string & name,
+                                            const bool project,
+                                            const libMesh::ParallelType type)
+{
+  if (_nl)
+    return &_nl->addVector(name, project, type);
+  else
+    return nullptr;
+}
diff --git a/framework/src/timeintegrators/TimeIntegrator.C b/framework/src/timeintegrators/TimeIntegrator.C
index 75ace47ad6ac..64d6999e6edf 100644
--- a/framework/src/timeintegrators/TimeIntegrator.C
+++ b/framework/src/timeintegrators/TimeIntegrator.C
@@ -30,13 +30,11 @@ TimeIntegrator::validParams()
 TimeIntegrator::TimeIntegrator(const InputParameters & parameters)
   : MooseObject(parameters),
     Restartable(this, "TimeIntegrators"),
+    NonlinearTimeIntegratorInterface(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base"),
+                                     *getCheckedPointerParam<SystemBase *>("_sys")),
+    LinearTimeIntegratorInterface(*getCheckedPointerParam<SystemBase *>("_sys")),
     _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
     _sys(*getCheckedPointerParam<SystemBase *>("_sys")),
-    _nl(_fe_problem.getNonlinearSystemBase(
-        dynamic_cast<NonlinearSystemBase *>(&_sys) ? _sys.number() : 0)),
-    _nonlinear_implicit_system(dynamic_cast<NonlinearImplicitSystem *>(&_sys.system())),
-    _Re_time(_nl.getResidualTimeVector()),
-    _Re_non_time(_nl.getResidualNonTimeVector()),
     _du_dot_du(_sys.duDotDus()),
     _solution(_sys.currentSolution()),
     _solution_old(_sys.solutionState(1)),
@@ -50,8 +48,6 @@ TimeIntegrator::TimeIntegrator(const InputParameters & parameters)
     _n_nonlinear_iterations(0),
     _n_linear_iterations(0),
     _is_lumped(false),
-    _u_dot_factor_tag(_fe_problem.addVectorTag("u_dot_factor", Moose::VECTOR_TAG_SOLUTION)),
-    _u_dotdot_factor_tag(_fe_problem.addVectorTag("u_dotdot_factor", Moose::VECTOR_TAG_SOLUTION)),
     _var_restriction(declareRestartableData<bool>(
         "var_restriction", !getParam<std::vector<VariableName>>("variables").empty())),
     _local_indices(declareRestartableData<std::vector<dof_id_type>>("local_indices")),
diff --git a/modules/solid_mechanics/src/timeintegrators/DirectCentralDifference.C b/modules/solid_mechanics/src/timeintegrators/DirectCentralDifference.C
index 228eed5398ed..40cbff69613f 100644
--- a/modules/solid_mechanics/src/timeintegrators/DirectCentralDifference.C
+++ b/modules/solid_mechanics/src/timeintegrators/DirectCentralDifference.C
@@ -90,35 +90,35 @@ DirectCentralDifference::solve()
   _nonlinear_implicit_system->update();
 
   // Calculating the lumped mass matrix for use in residual calculation
-  mass_matrix.vector_mult(_mass_matrix_diag, *_ones);
+  mass_matrix.vector_mult(*_mass_matrix_diag, *_ones);
 
   // Compute the residual
-  _explicit_residual.zero();
+  _explicit_residual->zero();
   _fe_problem.computeResidual(
-      *_nonlinear_implicit_system->current_local_solution, _explicit_residual, _nl.number());
+      *_nonlinear_implicit_system->current_local_solution, *_explicit_residual, _nl->number());
 
   // Move the residual to the RHS
-  _explicit_residual *= -1.0;
+  *_explicit_residual *= -1.0;
 
   // Perform the linear solve
   bool converged = performExplicitSolve(mass_matrix);
-  _nl.overwriteNodeFace(*_nonlinear_implicit_system->solution);
+  _nl->overwriteNodeFace(*_nonlinear_implicit_system->solution);
 
   // Update the solution
-  *_nonlinear_implicit_system->solution = _nl.solutionOld();
-  *_nonlinear_implicit_system->solution += _solution_update;
+  *_nonlinear_implicit_system->solution = _nl->solutionOld();
+  *_nonlinear_implicit_system->solution += *_solution_update;
 
   _nonlinear_implicit_system->update();
 
-  _nl.setSolution(*_nonlinear_implicit_system->current_local_solution);
+  _nl->setSolution(*_nonlinear_implicit_system->current_local_solution);
   _nonlinear_implicit_system->nonlinear_solver->converged = converged;
 }
 
 void
 DirectCentralDifference::postResidual(NumericVector<Number> & residual)
 {
-  residual += _Re_time;
-  residual += _Re_non_time;
+  residual += *_Re_time;
+  residual += *_Re_non_time;
   residual.close();
 
   // Reset time to the time at which to evaluate nodal BCs, which comes next
@@ -131,10 +131,10 @@ DirectCentralDifference::performExplicitSolve(SparseMatrix<Number> &)
   bool converged = false;
 
   // "Invert" the diagonal mass matrix
-  _mass_matrix_diag.reciprocal();
+  _mass_matrix_diag->reciprocal();
   // Calculate acceleration
   auto & accel = *_sys.solutionUDotDot();
-  accel.pointwise_mult(_mass_matrix_diag, _explicit_residual);
+  accel.pointwise_mult(*_mass_matrix_diag, *_explicit_residual);
 
   // Scaling the acceleration
   auto accel_scaled = accel.clone();
@@ -147,11 +147,11 @@ DirectCentralDifference::performExplicitSolve(SparseMatrix<Number> &)
   vel = *old_vel->clone();
   vel += *accel_scaled;
 
-  _solution_update = vel;
-  _solution_update.scale(_dt);
+  *_solution_update = vel;
+  _solution_update->scale(_dt);
 
   // Check for convergence by seeing if there is a nan or inf
-  auto sum = _solution_update.sum();
+  auto sum = _solution_update->sum();
   converged = std::isfinite(sum);
 
   // The linear iteration count remains zero
diff --git a/test/tests/time_integrators/implicit-euler/gold/ie-linearfv_out.e b/test/tests/time_integrators/implicit-euler/gold/ie-linearfv_out.e
new file mode 100644
index 000000000000..d1b96f4eac95
Binary files /dev/null and b/test/tests/time_integrators/implicit-euler/gold/ie-linearfv_out.e differ
diff --git a/test/tests/time_integrators/implicit-euler/ie-linearfv.i b/test/tests/time_integrators/implicit-euler/ie-linearfv.i
new file mode 100644
index 000000000000..9e2d36b86ab0
--- /dev/null
+++ b/test/tests/time_integrators/implicit-euler/ie-linearfv.i
@@ -0,0 +1,97 @@
+###########################################################
+# This is a simple test with a time-dependent problem
+# demonstrating the use of the TimeIntegrator system.
+#
+# Testing a solution that is second order in space
+# and first order in time
+#
+# @Requirement F1.30
+###########################################################
+
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  xmin = -1
+  xmax = 1
+  ymin = -1
+  ymax = 1
+  nx = 10
+  ny = 10
+[]
+
+[Problem]
+  linear_sys_names = 'u_sys'
+[]
+
+[Variables]
+  [u]
+    type = MooseLinearVariableFVReal
+    solver_sys = 'u_sys'
+    initial_condition = 0.0
+  []
+[]
+
+[Functions]
+  [forcing_fn]
+    type = ParsedFunction
+    expression = ((x*x)+(y*y))-(4*t)
+  []
+  [exact_fn]
+    type = ParsedFunction
+    expression = t*((x*x)+(y*y))
+  []
+[]
+
+[LinearFVKernels]
+  [ie]
+    type = LinearFVTimeDerivative
+    variable = u
+  []
+  [diff]
+    type = LinearFVDiffusion
+    variable = u
+  []
+  [source]
+    type = LinearFVSource
+    variable = u
+    source_density = forcing_fn
+  []
+[]
+
+[LinearFVBCs]
+  [all]
+    type = LinearFVAdvectionDiffusionFunctorDirichletBC
+    variable = u
+    boundary = '0 1 2 3'
+    functor = exact_fn
+  []
+[]
+
+[Postprocessors]
+  [l2_err]
+    type = ElementL2Error
+    variable = u
+    function = exact_fn
+  []
+[]
+
+[Executioner]
+  type = Transient
+
+  system_names = u_sys
+  l_tol = 1e-10
+
+  petsc_options_iname = '-pc_type -pc_hypre_type'
+  petsc_options_value = 'hypre boomeramg'
+
+  # Test of the TimeIntegrator System
+  scheme = 'implicit-euler'
+
+  start_time = 0.0
+  num_steps = 5
+  dt = 0.25
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/test/tests/time_integrators/implicit-euler/tests b/test/tests/time_integrators/implicit-euler/tests
index 2ead2c8c068c..d897a7af61a9 100644
--- a/test/tests/time_integrators/implicit-euler/tests
+++ b/test/tests/time_integrators/implicit-euler/tests
@@ -1,9 +1,9 @@
 [Tests]
   design = 'ImplicitEuler.md'
-  issues = '#1953'
+  issues = '#1953 #29149'
 
   [group]
-    requirement = "The system shall support the use of an implicit Euler solver with"
+    requirement = "The system shall support the use of an implicit Euler solver "
 
     [test]
       type = 'Exodiff'
@@ -11,7 +11,7 @@
       exodiff = 'ie_out.e'
       use_old_floor = true
       abs_zero = 1e-9
-      detail = "with and"
+      detail = "without and"
     []
 
     [adapt]
@@ -19,10 +19,19 @@
       input = 'ie_adapt.i'
       exodiff = 'ie_adapt_out.e-s005'
       group = 'adaptive'
-      detail = "without mesh adaptivity."
+      detail = "with mesh adaptivity."
     []
   []
 
+  [linearfv]
+    type = 'Exodiff'
+    input = 'ie-linearfv.i'
+    exodiff = 'ie-linearfv_out.e'
+    abs_zero = 1e-9
+    requirement = "The system shall be able to use implicit euler time integration with linear FV spatial discretization."
+    max_threads = 1
+  []
+
   [monomials]
     requirement = "The system shall support the use of an implicit Euler solver with discontinuous "
                   "(first-order Monomial) shape functions."