diff --git a/src/loMach.cpp b/src/loMach.cpp
index 337efd20e..9f95926f7 100644
--- a/src/loMach.cpp
+++ b/src/loMach.cpp
@@ -114,6 +114,7 @@ void LoMachSolver::initialize() {
     iter = 0;
     temporal_coeff_.nStep = iter;
   }
+  temporal_coeff_.order = max_bdf_order;
 
   //-----------------------------------------------------
   // 1) Prepare the mesh
diff --git a/src/loMach.hpp b/src/loMach.hpp
index dc4efc38c..ba5a9066a 100644
--- a/src/loMach.hpp
+++ b/src/loMach.hpp
@@ -74,6 +74,9 @@ class Tps;
 #include "turb_model_base.hpp"
 
 struct temporalSchemeCoefficients {
+  // Max order
+  int order;
+
   // Current time
   double time;
 
diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp
index d389f78cf..54e8b2234 100644
--- a/src/reactingFlow.cpp
+++ b/src/reactingFlow.cpp
@@ -421,6 +421,23 @@ ReactingFlow::ReactingFlow(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, tem
   chemistry_ = new Chemistry(mixture_, chemistryInput_);
 
   tpsP_->getInput("loMach/reactingFlow/ic", ic_string_, std::string(""));
+  tpsP_->getInput("loMach/reactingFlow/sub-steps", nSub_, 1);
+
+  // Check time marching order.  Operator split (i.e., nSub_ > 1) not supported for order > 1.
+  if ((nSub_ > 1) && (time_coeff_.order > 1)) {
+    if (rank0_) {
+      std::cout << "ERROR: BDF order > 1 not supported with operator split." << std::endl;
+      std::cout << "       Either set loMach/reactingFlow/sub-steps = 1 or time/bdfOrder = 1" << std::endl;
+    }
+    assert(false);
+    exit(1);
+  }
+
+  // By default, do not use operator split scheme.
+  operator_split_ = false;
+  if (nSub_ > 1) {
+    operator_split_ = true;
+  }
 }
 
 ReactingFlow::~ReactingFlow() {
@@ -566,6 +583,12 @@ void ReactingFlow::initializeSelf() {
   NYnm1_ = 0.0;
   NYnm2_ = 0.0;
 
+  // for time-splitting
+  spec_buffer_.SetSize(yDofInt_);
+  temp_buffer_.SetSize(sDofInt_);
+  YnStar_.SetSize(yDofInt_);
+  TnStar_.SetSize(sDofInt_);
+
   // number density
   Xn_.SetSize(yDofInt_);
   Xn_ = 0.0;
@@ -1235,10 +1258,25 @@ void ReactingFlow::step() {
   updateMixture();
   extrapolateState();
   updateBC(0);  // NB: can't ramp right now
-  updateThermoP();
+
+  if (!operator_split_) {
+    updateThermoP();
+  }
+
   updateDensity(1.0);
   updateDiffusivity();
-  speciesProduction();
+
+  // PART I: Form and solve implicit systems for species and temperature.
+  //
+  // If NOT using operator splitting, these systems include all terms,
+  // with the reaction source terms treated explicitly in time.  When
+  // using operator splitting, this solve includes the
+  // advection-diffusion terms but the reaction source terms are
+  // handled separately with substepping below.
+
+  if (!operator_split_) {
+    speciesProduction();
+  }
 
   // advance species, last slot is from calculated sum of others
   for (int iSpecies = 0; iSpecies < nSpecies_ - 1; iSpecies++) {
@@ -1248,9 +1286,10 @@ void ReactingFlow::step() {
   YnFull_gf_.SetFromTrueDofs(Yn_next_);
 
   // ho_i*w_i
-  heatOfFormation();
+  if (!operator_split_) {
+    heatOfFormation();
+  }
 
-  // TODO(swh): this is a bad name
   crossDiffusion();
 
   // advance temperature
@@ -1268,11 +1307,55 @@ void ReactingFlow::step() {
     Tn_next_gf_.GetTrueDofs(Tn_next_);
   }
 
-  // explicit filter species (is this necessary?)
+  // explicit filter species
+  // Add if found to be necessary
+
+  if (operator_split_) {
+    /// PART II: time-splitting of reaction
+
+    // Save Yn_ and Tn_ because substep routines overwrite these as
+    // with the {n}+iSub substates, which is convenient b/c helper
+    // functions (like speciesProduction) use these Vectors
+    spec_buffer_.Set(1.0, Yn_);
+    temp_buffer_.Set(1.0, Tn_);
 
-  // prepare for external use
+    // Evaluate (Yn_next_ - Yn_)/nSub and store in YnStar_, and analog
+    // for TnStar_.
+    substepState();
+
+    for (int iSub = 0; iSub < nSub_; iSub++) {
+      // update wdot quantities at full substep in Yn/Tn state
+      updateMixture();
+      updateThermoP();
+      updateDensity(0.0);
+      speciesProduction();
+      heatOfFormation();
+
+      // advance over substep
+      for (int iSpecies = 0; iSpecies < nSpecies_ - 1; iSpecies++) {
+        speciesSubstep(iSpecies, iSub);
+      }
+      speciesLastSubstep();
+      Yn_gf_.SetFromTrueDofs(Yn_);
+
+      temperatureSubstep(iSub);
+      Tn_gf_.SetFromTrueDofs(Tn_);
+    }
+
+    // set register to correct full step fields
+    Yn_next_ = Yn_;
+    Tn_next_ = Tn_;
+
+    YnFull_gf_.SetFromTrueDofs(Yn_next_);
+    Tn_next_gf_.SetFromTrueDofs(Tn_next_);
+
+    Yn_ = spec_buffer_;
+    Tn_ = temp_buffer_;
+    Tn_gf_.SetFromTrueDofs(Tn_);
+  }
+
+  /// PART III: prepare for external use
   updateDensity(1.0);
-  // computeQt();
   computeQtTO();
 
   UpdateTimestepHistory(dt_);
@@ -1298,12 +1381,17 @@ void ReactingFlow::temperatureStep() {
   MsRhoCp_form_->FormSystemMatrix(empty, MsRhoCp_);
   MsRhoCp_->AddMult(tmpR0_, resT_, -1.0);
 
-  // dPo/dt
-  tmpR0_ = dtP_;  // with rho*Cp on LHS, no Cp on this term
-  Ms_->AddMult(tmpR0_, resT_);
+  // If NOT using operator splitting, include the unsteady pressure
+  // term and the heat of formation contribution on the rhs.  If using
+  // operator splitting, these are handled in temperatureSubstep.
+  if (!operator_split_) {
+    // dPo/dt
+    tmpR0_ = dtP_;  // with rho*Cp on LHS, no Cp on this term
+    Ms_->AddMult(tmpR0_, resT_);
 
-  // heat of formation
-  Ms_->AddMult(hw_, resT_);
+    // heat of formation
+    Ms_->AddMult(hw_, resT_);
+  }
 
   // species-temp diffusion term, already in int-weak form
   resT_.Add(1.0, crossDiff_);
@@ -1351,6 +1439,31 @@ void ReactingFlow::temperatureStep() {
   Tn_next_gf_.GetTrueDofs(Tn_next_);
 }
 
+void ReactingFlow::temperatureSubstep(int iSub) {
+  // substep dt
+  double dtSub = dt_ / (double)nSub_;
+
+  CpMix_gf_.GetTrueDofs(CpMix_);
+
+  // heat of formation term
+  tmpR0_.Set(1.0, hw_);
+
+  // pressure
+  tmpR0_ += dtP_;
+
+  // contribution to temperature update is dt * rhs / (rho * Cp)
+  tmpR0_ /= rn_;
+  tmpR0_ /= CpMix_;
+  tmpR0_ *= dtSub;
+
+  // TnStar has star state at substep here
+  tmpR0_.Add(1.0, TnStar_);
+  tmpR0_.Add(1.0, Tn_);
+
+  // Tn now has full state at substep
+  Tn_.Set(1.0, tmpR0_);
+}
+
 void ReactingFlow::speciesLastStep() {
   tmpR0_ = 0.0;
   double *dataYn = Yn_next_.HostReadWrite();
@@ -1365,6 +1478,20 @@ void ReactingFlow::speciesLastStep() {
   }
 }
 
+void ReactingFlow::speciesLastSubstep() {
+  tmpR0_ = 0.0;
+  double *dataYn = Yn_.HostReadWrite();
+  double *dataYsum = tmpR0_.HostReadWrite();
+  for (int n = 0; n < nSpecies_ - 1; n++) {
+    for (int i = 0; i < sDofInt_; i++) {
+      dataYsum[i] += dataYn[i + n * sDofInt_];
+    }
+  }
+  for (int i = 0; i < sDofInt_; i++) {
+    dataYn[i + (nSpecies_ - 1) * sDofInt_] = 1.0 - dataYsum[i];
+  }
+}
+
 void ReactingFlow::speciesStep(int iSpec) {
   Array<int> empty;
 
@@ -1392,9 +1519,13 @@ void ReactingFlow::speciesStep(int iSpec) {
   MsRho_form_->FormSystemMatrix(empty, MsRho_);
   MsRho_->AddMult(tmpR0_, resY_, -1.0);
 
-  // production of iSpec
-  setScalarFromVector(prodY_, iSpec, &tmpR0_);
-  Ms_->AddMult(tmpR0_, resY_);
+  // If NOT using operator splitting, include the reaction source term
+  // on the rhs.  If using operator splitting, this is handled by
+  // speciesSubstep.
+  if (!operator_split_) {
+    setScalarFromVector(prodY_, iSpec, &tmpR0_);
+    Ms_->AddMult(tmpR0_, resY_);
+  }
 
   // Add natural boundary terms here later
   // NB: adiabatic natural BC is handled, but don't have ability to impose non-zero heat flux yet
@@ -1445,6 +1576,27 @@ void ReactingFlow::speciesStep(int iSpec) {
   setVectorFromScalar(tmpR0_, iSpec, &Yn_next_);
 }
 
+// Y{n + (substep+1)} = dt * wdot(Y{n + (substep)} + Y{n + (substep+1)}*
+void ReactingFlow::speciesSubstep(int iSpec, int iSub) {
+  // substep dt
+  double dtSub = dt_ / (double)nSub_;
+
+  // production of iSpec
+  setScalarFromVector(prodY_, iSpec, &tmpR0_);
+  tmpR0_ /= rn_;
+  tmpR0_ *= dtSub;
+
+  // YnStar has star state at substep here
+  setScalarFromVector(YnStar_, iSpec, &tmpR0a_);
+  tmpR0_.Add(1.0, tmpR0a_);
+
+  setScalarFromVector(Yn_, iSpec, &tmpR0a_);
+  tmpR0_.Add(1.0, tmpR0a_);
+
+  // Yn now has full state at substep
+  setVectorFromScalar(tmpR0_, iSpec, &Yn_);
+}
+
 void ReactingFlow::speciesProduction() {
   const double *dataT = Tn_.HostRead();
   const double *dataRho = rn_.HostRead();
@@ -1657,6 +1809,19 @@ void ReactingFlow::extrapolateState() {
   Yext_.Add(time_coeff_.ab3, Ynm2_);
 }
 
+// Interpolate star state to substeps. Linear for now, use higher order
+// if time-splitting order can be improved
+void ReactingFlow::substepState() {
+  // delta over substep of star state
+  double wt = 1.0 / (double)nSub_;
+
+  TnStar_.Set(wt, Tn_next_);
+  TnStar_.Add(-wt, temp_buffer_);
+
+  YnStar_.Set(wt, Yn_next_);
+  YnStar_.Add(-wt, spec_buffer_);
+}
+
 void ReactingFlow::updateMixture() {
   {
     double *dataY = Yn_gf_.HostReadWrite();
@@ -1765,7 +1930,11 @@ void ReactingFlow::updateThermoP() {
     thermo_pressure_ = system_mass_ / allMass * Pnm1_;
     dtP_ = time_coeff_.bd0 * thermo_pressure_ + time_coeff_.bd1 * Pnm1_ + time_coeff_.bd2 * Pnm2_ +
            time_coeff_.bd3 * Pnm3_;
-    dtP_ *= 1.0 / dt_;
+
+    // For operator split, dtP term is in split substeps, so need to
+    // divide dt by nSub.  For 'unified' approach (not operator
+    // split), nSub = 1, so it makes no difference.
+    dtP_ *= (nSub_ / dt_);
   }
 }
 
diff --git a/src/reactingFlow.hpp b/src/reactingFlow.hpp
index 854794215..80e734ca9 100644
--- a/src/reactingFlow.hpp
+++ b/src/reactingFlow.hpp
@@ -282,6 +282,7 @@ class ReactingFlow : public ThermoChemModelBase {
   Vector specificHeatRatios_;
   Vector initialMassFraction_;
   Vector atomMW_;
+  Vector CpMix_;
 
   Vector Qt_;
   Vector rn_;
@@ -289,6 +290,12 @@ class ReactingFlow : public ThermoChemModelBase {
   Vector kappa_;
   Vector visc_;
 
+  // time-splitting
+  Vector YnStar_, spec_buffer_;
+  Vector TnStar_, temp_buffer_;
+  bool operator_split_ = false;
+  int nSub_;
+
   // Parameters and objects used in filter-based stabilization
   bool filter_temperature_ = false;
   int filter_cutoff_modes_ = 0;
@@ -331,6 +338,12 @@ class ReactingFlow : public ThermoChemModelBase {
   void heatOfFormation();
   void crossDiffusion();
 
+  // time-splitting
+  void substepState();
+  void speciesLastSubstep();
+  void speciesSubstep(int iSpec, int iSub);
+  void temperatureSubstep(int iSub);
+
   /// for creation of structs to interface with old plasma/chem stuff
   void identifySpeciesType(Array<ArgonSpcs> &speciesType);
   void identifyCollisionType(const Array<ArgonSpcs> &speciesType, ArgonColl *collisionIndex);
diff --git a/test/inputs/input.reactStiffRx.ini b/test/inputs/input.reactStiffRx.ini
new file mode 100644
index 000000000..b2120984c
--- /dev/null
+++ b/test/inputs/input.reactStiffRx.ini
@@ -0,0 +1,137 @@
+[solver]
+#type = flow
+type = loMach
+
+[loMach]
+mesh = meshes/periodic-square.mesh
+order = 2
+maxIters = 50
+outputFreq = 1
+timingFreq = 1
+fluid = 'user_defined'
+#viscosityMultiplier = 99000.
+#equation_system = 'navier-stokes'
+enablePressureForcing = false
+pressureGrad = '0.0 0.0 0.0'
+enableGravity = false
+gravity = '0.0 0.0 0.0'
+openSystem = False
+ambientPressure = 101300. # to match input.reaction.ini
+flow-solver = tomboulides
+thermo-solver = reacting-flow
+
+[loMach/tomboulides]
+linear-solver-rtol = 1.0e-12
+linear-solver-max-iter = 2000
+
+[loMach/reactingFlow]
+sub-steps = 20
+linear-solver-rtol = 1.0e-12
+linear-solver-max-iter = 2000
+initialTemperature = 294.075;
+# ic = 'step'
+
+[io]
+outdirBase = output-stiff
+
+[time]
+integrator = curlcurl
+dt_fixed = 2.0e-4
+bdfOrder = 1
+
+[spongeMultiplier]
+uniform = true
+uniformMult = 99000.
+
+[initialConditions]
+temperature = 294.0
+
+[boundaryConditions]
+numWalls = 0
+numInlets = 0
+numOutlets = 0
+
+#[periodicity]
+#enablePeriodic = True
+#periodicX = True
+#periodicY = True
+#periodicZ = True
+
+###########################################
+#        PLASMA MODELS
+###########################################
+[plasma_models]
+ambipolar = False
+two_temperature = False
+gas_model = perfect_mixture
+transport_model = constant
+chemistry_model = 'mass_action_law'
+
+[plasma_models/transport_model/constant]
+viscosity = 0.0
+bulk_viscosity = 0.0
+thermal_conductivity = 0.0
+electron_thermal_conductivity = 0.0
+diffusivity/species1 = 0.0
+diffusivity/species2 = 0.0
+diffusivity/species3 = 0.0
+momentum_transfer_frequency/species1 = 0.0
+momentum_transfer_frequency/species2 = 0.0
+momentum_transfer_frequency/species3 = 0.0
+
+[atoms]
+numAtoms = 2
+
+[atoms/atom1]
+name = 'Ar'
+mass = 2.896439e-2
+
+[atoms/atom2]
+name = 'E'
+mass = 1e-7
+
+[species]
+numSpecies = 3
+background_index = 3
+
+[species/species1]
+name = 'Ar.+1'
+composition = '{Ar : 1, E : -1}'
+#initialMassFraction = 0.
+initialMassFraction = 1.0e-12
+formation_energy = 10000
+perfect_mixture/constant_molar_cv = 2.49996
+
+[species/species2]
+name = 'E'
+composition = '{E : 1}'
+#initialMassFraction = 0.
+initialMassFraction = 1.0e-12
+initialElectronTemperature = 300.;
+formation_energy = 0.
+perfect_mixture/constant_molar_cv = 2.49996
+
+[species/species3]
+name = 'Ar'
+composition = '{Ar : 1}'
+initialMassFraction = 1.
+formation_energy = 0.
+perfect_mixture/constant_molar_cv = 2.49996
+
+[reactions]
+number_of_reactions = 1
+
+[reactions/reaction1]
+equation = 'Ar + E <=> Ar.+1 + 2 E'
+reaction_energy = 1.0e4
+reactant_stoichiometry = '0 0 1'
+product_stoichiometry = '1 1 0'
+model = arrhenius
+arrhenius/A = 1e-7
+arrhenius/b = 4.
+arrhenius/E = 0.
+
+detailed_balance = True
+equilibrium_constant/A = 1e-10
+equilibrium_constant/b = 4.
+equilibrium_constant/E = 0.
diff --git a/test/reactFlow-singleRx.test b/test/reactFlow-singleRx.test
index 99f32d8dd..d925494e5 100755
--- a/test/reactFlow-singleRx.test
+++ b/test/reactFlow-singleRx.test
@@ -3,12 +3,16 @@
 
 TEST="singleReaction"
 RUNFILE="inputs/input.reactSingleRx.ini"
+RUNFILE_STIFF="inputs/input.reactStiffRx.ini"
 EXE="../src/tps"
 RESTART="ref_solns/reactSingleRx/restart_output.sol.h5"
 
 setup() {
     SOLN_FILE=restart_output.sol.h5
     REF_FILE=ref_solns/reactSingleRx/restart_output.sol.h5
+
+    SOLN_FILE_STIFF=restart_output-stiff.sol.h5
+    REF_FILE_STIFF=ref_solns/reactSingleRx/restart_output-stiff.sol.h5
 }
 
 @test "[$TEST] check for input file $RUNFILE" {
@@ -19,7 +23,7 @@ setup() {
     rm -rf output/*
     rm -f $SOLN_FILE
     $EXE --runFile $RUNFILE
-    test -s $SOLN_FILE    
+    test -s $SOLN_FILE
 }
 
 @test "[$TEST] verify tps output with input -> $RUNFILE" {
@@ -27,3 +31,21 @@ setup() {
     test -s $REF_FILE
     h5diff -r --delta=1e-11 $SOLN_FILE $REF_FILE /species/speciesAll
 }
+
+@test "[$TEST] check for input file $RUNFILE_STIFF" {
+    test -s $RUNFILE
+}
+
+@test "[$TEST] run tps with input -> $RUNFILE_STIFF" {
+    rm -rf output-stiff
+    rm -f $SOLN_FILE_STIFF
+    $EXE --runFile $RUNFILE_STIFF
+    test -s $SOLN_FILE_STIFF
+}
+
+@test "[$TEST] verify tps output with input -> $RUNFILE_STIFF" {
+    test -s $SOLN_FILE_STIFF
+    test -s $REF_FILE_STIFF
+    h5diff -r --relative=1e-12 $SOLN_FILE_STIFF $REF_FILE_STIFF /species/speciesAll
+    h5diff -r --relative=1e-12 $SOLN_FILE_STIFF $REF_FILE_STIFF /temperature/temperature
+}
diff --git a/test/ref_solns/.gitattributes b/test/ref_solns/.gitattributes
index d87ce5fe7..0aaa8b015 100644
--- a/test/ref_solns/.gitattributes
+++ b/test/ref_solns/.gitattributes
@@ -6,3 +6,4 @@ ref_solns/heatedBox/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -text
 taylor-couette/restart_output-tc.Re100.h5 filter=lfs diff=lfs merge=lfs -text
 pipe/restart_output-pipe-lam.sol.h5 filter=lfs diff=lfs merge=lfs -text
 pipe/restart_output-pipe-arans.sol.h5 filter=lfs diff=lfs merge=lfs -text
+reactSingleRx/restart_output-stiff.sol.h5 filter=lfs diff=lfs merge=lfs -text
diff --git a/test/ref_solns/reactSingleRx/restart_output-stiff.sol.h5 b/test/ref_solns/reactSingleRx/restart_output-stiff.sol.h5
new file mode 100644
index 000000000..3f1cb2f5d
--- /dev/null
+++ b/test/ref_solns/reactSingleRx/restart_output-stiff.sol.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fc10a16abb8c1e2cf0725948c85f083022bc3d09e24c48139735c2a2be0d58cd
+size 7896