diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h
index 9509d842b..5ebfbdcd9 100644
--- a/include/cantera/oneD/StFlow.h
+++ b/include/cantera/oneD/StFlow.h
@@ -288,6 +288,16 @@ public:
     virtual void evalContinuity(size_t j, doublereal* x, doublereal* r,
                                 integer* diag, doublereal rdt) = 0;
 
+    //! Index of the species on the left boundary with the largest mass fraction
+    size_t leftExcessSpecies() const {
+        return m_kExcessLeft;
+    }
+
+    //! Index of the species on the right boundary with the largest mass fraction
+    size_t rightExcessSpecies() const {
+        return m_kExcessRight;
+    }
+
 protected:
     doublereal component(const doublereal* x, size_t i, size_t j) const {
         return x[index(i,j)];
@@ -487,6 +497,12 @@ protected:
     vector_fp m_zfix;
     vector_fp m_tfix;
 
+    //! Index of species with a large mass fraction at each boundary, for which
+    //! the mass fraction may be calculated as 1 minus the sum of the other mass
+    //! fractions
+    size_t m_kExcessLeft;
+    size_t m_kExcessRight;
+
     bool m_dovisc;
 
     //! Update the transport properties at grid points in the range from `j0`
diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp
index 77209488d..9a9ee138f 100644
--- a/src/oneD/StFlow.cpp
+++ b/src/oneD/StFlow.cpp
@@ -24,7 +24,9 @@ StFlow::StFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
     m_epsilon_right(0.0),
     m_do_soret(false),
     m_transport_option(-1),
-    m_do_radiation(false)
+    m_do_radiation(false),
+    m_kExcessLeft(0),
+    m_kExcessRight(0)
 {
     m_type = cFlowType;
     m_points = points;
@@ -262,9 +264,14 @@ void StFlow::eval(size_t jg, doublereal* xg,
     // ------------ update properties ------------
 
     updateThermo(x, j0, j1);
-    // update transport properties only if a Jacobian is not being evaluated
     if (jg == npos) {
+        // update transport properties only if a Jacobian is not being evaluated
         updateTransport(x, j0, j1);
+
+        double* Yleft = x + index(c_offset_Y, jmin);
+        m_kExcessLeft = distance(Yleft, max_element(Yleft, Yleft + m_nsp));
+        double* Yright = x + index(c_offset_Y, jmax);
+        m_kExcessRight = distance(Yright, max_element(Yright, Yright + m_nsp));
     }
 
     // update the species diffusive mass fluxes whether or not a
@@ -373,7 +380,7 @@ void StFlow::eval(size_t jg, doublereal* xg,
                 rsd[index(c_offset_Y + k, 0)] =
                     -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0));
             }
-            rsd[index(c_offset_Y, 0)] = 1.0 - sum;
+            rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum;
         } else if (j == m_points - 1) {
             evalRightBoundary(x, rsd, diag, rdt);
         } else { // interior points
@@ -871,8 +878,8 @@ void AxiStagnFlow::evalRightBoundary(doublereal* x, doublereal* rsd,
         sum += Y(x,k,j);
         rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j);
     }
-    rsd[index(4,j)] = 1.0 - sum;
-    diag[index(4,j)] = 0;
+    rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum;
+    diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0;
 }
 
 void AxiStagnFlow::evalContinuity(size_t j, doublereal* x, doublereal* rsd,
@@ -924,8 +931,8 @@ void FreeFlame::evalRightBoundary(doublereal* x, doublereal* rsd,
         sum += Y(x,k,j);
         rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j);
     }
-    rsd[index(4,j)] = 1.0 - sum;
-    diag[index(4,j)] = 0;
+    rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum;
+    diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0;
 }
 
 void FreeFlame::evalContinuity(size_t j, doublereal* x, doublereal* rsd,
diff --git a/src/oneD/boundaries1D.cpp b/src/oneD/boundaries1D.cpp
index 42eaf3af5..ae37fcb42 100644
--- a/src/oneD/boundaries1D.cpp
+++ b/src/oneD/boundaries1D.cpp
@@ -185,8 +185,10 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg,
         rb[3] += x[0];
 
         // add the convective term to the species residual equations
-        for (size_t k = 1; k < m_nsp; k++) {
-            rb[4+k] += x[0]*m_yin[k];
+        for (size_t k = 0; k < m_nsp; k++) {
+            if (k != m_flow_right->leftExcessSpecies()) {
+                rb[4+k] += x[0]*m_yin[k];
+            }
         }
 
         // if the flow is a freely-propagating flame, mdot is not specified.
@@ -204,8 +206,10 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg,
         rb[1] -= m_V0;
         rb[2] -= x[1]; // T
         rb[0] += x[0]; // u
-        for (size_t k = 1; k < m_nsp; k++) {
-            rb[4+k] += x[0]*m_yin[k];
+        for (size_t k = 0; k < m_nsp; k++) {
+            if (k != m_flow_left->rightExcessSpecies()) {
+                rb[4+k] += x[0]*m_yin[k];
+            }
         }
     }
 }
@@ -444,9 +448,12 @@ void Outlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg,
         }
 
         rb[2] = xb[2] - xb[2 - nc]; // zero T gradient
-        for (k = 5; k < nc; k++) {
-            rb[k] = xb[k] - xb[k - nc]; // zero mass fraction gradient
-            db[k] = 0;
+        size_t kSkip = 4 + m_flow_left->rightExcessSpecies();
+        for (k = 4; k < nc; k++) {
+            if (k != kSkip) {
+                rb[k] = xb[k] - xb[k - nc]; // zero mass fraction gradient
+                db[k] = 0;
+            }
         }
     }
 }
@@ -573,9 +580,12 @@ void OutletRes1D::eval(size_t jg, doublereal* xg, doublereal* rg,
             rb[0] = xb[3]; // zero Lambda
         }
         rb[2] = xb[2] - m_temp; // zero dT/dz
-        for (k = 5; k < nc; k++) {
-            rb[k] = xb[k] - m_yres[k-4]; // fixed Y
-            db[k] = 0;
+        size_t kSkip = m_flow_left->rightExcessSpecies();
+        for (k = 4; k < nc; k++) {
+            if (k != kSkip) {
+                rb[k] = xb[k] - m_yres[k-4]; // fixed Y
+                db[k] = 0;
+            }
         }
     }
 }
@@ -795,8 +805,11 @@ void ReactingSurf1D::eval(size_t jg, doublereal* xg, doublereal* rg,
         rb =r - nc;
         xb = x - nc;
         rb[2] = xb[2] - x[0]; // specified T
-        for (size_t nl = 1; nl < m_left_nsp; nl++) {
-            rb[4+nl] += m_work[nl]*mwleft[nl];
+        size_t nSkip = m_flow_left->rightExcessSpecies();
+        for (size_t nl = 0; nl < m_left_nsp; nl++) {
+            if (nl != nSkip) {
+                rb[4+nl] += m_work[nl]*mwleft[nl];
+            }
         }
     }
 }