diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h
index d1e73030a..0d19460d9 100644
--- a/include/cantera/zeroD/ReactorBase.h
+++ b/include/cantera/zeroD/ReactorBase.h
@@ -57,7 +57,6 @@ public:
      */
     void setInitialVolume(doublereal vol) {
         m_vol = vol;
-        m_vol0 = vol;
     }
 
     /**
@@ -210,10 +209,8 @@ protected:
     size_t m_nsp;
 
     thermo_t*  m_thermo;
-    doublereal m_vol, m_vol0;
+    doublereal m_vol;
     bool m_init;
-    size_t m_nInlets, m_nOutlets;
-    bool m_open;
     doublereal m_enthalpy;
     doublereal m_intEnergy;
     doublereal m_pressure;
@@ -221,9 +218,7 @@ protected:
     std::vector<FlowDevice*> m_inlet, m_outlet;
     std::vector<Wall*> m_wall;
     vector_int m_lr;
-    size_t m_nwalls;
     std::string m_name;
-    double m_rho0;
 
     //! The ReactorNet that this reactor is part of
     ReactorNet* m_net;
diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp
index 9e427b284..142ed0a81 100644
--- a/src/zeroD/ConstPressureReactor.cpp
+++ b/src/zeroD/ConstPressureReactor.cpp
@@ -94,27 +94,23 @@ void ConstPressureReactor::evalEqs(doublereal time, doublereal* y,
     // external heat transfer
     double dHdt = - m_Q;
 
-    // add terms for open system
-    if (m_open) {
-        double enthalpy = m_thermo->enthalpy_mass();
-        // outlets
-        for (size_t i = 0; i < m_nOutlets; i++) {
-            double mdot_out = m_outlet[i]->massFlowRate(time); // mass flow out of system
-            dmdt -= mdot_out;
-            dHdt -= mdot_out * enthalpy;
-        }
+    // add terms for outlets
+    for (size_t i = 0; i < m_outlet.size(); i++) {
+        double mdot_out = m_outlet[i]->massFlowRate(time); // mass flow out of system
+        dmdt -= mdot_out;
+        dHdt -= mdot_out * m_enthalpy;
+    }
 
-        // inlets
-        for (size_t i = 0; i < m_nInlets; i++) {
-            double mdot_in = m_inlet[i]->massFlowRate(time);
-            dmdt += mdot_in; // mass flow into system
-            for (size_t n = 0; n < m_nsp; n++) {
-                double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
-                // flow of species into system and dilution by other species
-                dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
-            }
-            dHdt += mdot_in * m_inlet[i]->enthalpy_mass();
+    // add terms for inlets
+    for (size_t i = 0; i < m_inlet.size(); i++) {
+        double mdot_in = m_inlet[i]->massFlowRate(time);
+        dmdt += mdot_in; // mass flow into system
+        for (size_t n = 0; n < m_nsp; n++) {
+            double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
+            // flow of species into system and dilution by other species
+            dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
         }
+        dHdt += mdot_in * m_inlet[i]->enthalpy_mass();
     }
 
     ydot[0] = dmdt;
diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp
index 527f717ff..ee671af2c 100644
--- a/src/zeroD/IdealGasConstPressureReactor.cpp
+++ b/src/zeroD/IdealGasConstPressureReactor.cpp
@@ -106,24 +106,21 @@ void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y,
         dYdt[n] -= Y[n] * mdot_surf / m_mass;
     }
 
-    // add terms for open system
-    if (m_open) {
-        // outlets
-        for (size_t i = 0; i < m_nOutlets; i++) {
-            dmdt -= m_outlet[i]->massFlowRate(time); // mass flow out of system
-        }
+    // add terms for outlets
+    for (size_t i = 0; i < m_outlet.size(); i++) {
+        dmdt -= m_outlet[i]->massFlowRate(time); // mass flow out of system
+    }
 
-        // inlets
-        for (size_t i = 0; i < m_nInlets; i++) {
-            double mdot_in = m_inlet[i]->massFlowRate(time);
-            dmdt += mdot_in; // mass flow into system
-            mcpdTdt += m_inlet[i]->enthalpy_mass() * mdot_in;
-            for (size_t n = 0; n < m_nsp; n++) {
-                double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
-                // flow of species into system and dilution by other species
-                dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
-                mcpdTdt -= m_hk[n] / mw[n] * mdot_spec;
-            }
+    // add terms for inlets
+    for (size_t i = 0; i < m_inlet.size(); i++) {
+        double mdot_in = m_inlet[i]->massFlowRate(time);
+        dmdt += mdot_in; // mass flow into system
+        mcpdTdt += m_inlet[i]->enthalpy_mass() * mdot_in;
+        for (size_t n = 0; n < m_nsp; n++) {
+            double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
+            // flow of species into system and dilution by other species
+            dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
+            mcpdTdt -= m_hk[n] / mw[n] * mdot_spec;
         }
     }
 
diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp
index 50a79f27f..9ee21021b 100644
--- a/src/zeroD/IdealGasReactor.cpp
+++ b/src/zeroD/IdealGasReactor.cpp
@@ -115,29 +115,26 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
         dYdt[n] -= Y[n] * mdot_surf / m_mass;
     }
 
-    // add terms for open system
-    if (m_open) {
-        // outlets
-        for (size_t i = 0; i < m_nOutlets; i++) {
-            double mdot_out = m_outlet[i]->massFlowRate(time);
-            dmdt -= mdot_out; // mass flow out of system
-            mcvdTdt -= mdot_out * m_pressure * m_vol / m_mass; // flow work
-        }
+    // add terms for outlets
+    for (size_t i = 0; i < m_outlet.size(); i++) {
+        double mdot_out = m_outlet[i]->massFlowRate(time);
+        dmdt -= mdot_out; // mass flow out of system
+        mcvdTdt -= mdot_out * m_pressure * m_vol / m_mass; // flow work
+    }
 
-        // inlets
-        for (size_t i = 0; i < m_nInlets; i++) {
-            double mdot_in = m_inlet[i]->massFlowRate(time);
-            dmdt += mdot_in; // mass flow into system
-            mcvdTdt += m_inlet[i]->enthalpy_mass() * mdot_in;
-            for (size_t n = 0; n < m_nsp; n++) {
-                double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
-                // flow of species into system and dilution by other species
-                dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
+    // add terms for inlets
+    for (size_t i = 0; i < m_inlet.size(); i++) {
+        double mdot_in = m_inlet[i]->massFlowRate(time);
+        dmdt += mdot_in; // mass flow into system
+        mcvdTdt += m_inlet[i]->enthalpy_mass() * mdot_in;
+        for (size_t n = 0; n < m_nsp; n++) {
+            double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
+            // flow of species into system and dilution by other species
+            dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
 
-                // In combintion with h_in*mdot_in, flow work plus thermal
-                // energy carried with the species
-                mcvdTdt -= m_uk[n] / mw[n] * mdot_spec;
-            }
+            // In combintion with h_in*mdot_in, flow work plus thermal
+            // energy carried with the species
+            mcvdTdt -= m_uk[n] / mw[n] * mdot_spec;
         }
     }
 
diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp
index 3af7ffc4f..20939046b 100644
--- a/src/zeroD/Reactor.cpp
+++ b/src/zeroD/Reactor.cpp
@@ -56,7 +56,7 @@ void Reactor::getInitialConditions(double t0, size_t leny, double* y)
 void Reactor::getSurfaceInitialConditions(double* y)
 {
     size_t loc = 0;
-    for (size_t m = 0; m < m_nwalls; m++) {
+    for (size_t m = 0; m < m_wall.size(); m++) {
         SurfPhase* surf = m_wall[m]->surface(m_lr[m]);
         if (surf) {
             m_wall[m]->getCoverages(m_lr[m], y + loc);
@@ -71,7 +71,7 @@ void Reactor::initialize(doublereal t0)
     m_sdot.resize(m_nsp, 0.0);
     m_wdot.resize(m_nsp, 0.0);
     m_nv = m_nsp + 3;
-    for (size_t w = 0; w < m_nwalls; w++)
+    for (size_t w = 0; w < m_wall.size(); w++)
         if (m_wall[w]->surface(m_lr[w])) {
             m_nv += m_wall[w]->surface(m_lr[w])->nSpecies();
         }
@@ -81,7 +81,7 @@ void Reactor::initialize(doublereal t0)
     m_intEnergy = m_thermo->intEnergy_mass();
 
     size_t nt = 0, maxnt = 0;
-    for (size_t m = 0; m < m_nwalls; m++) {
+    for (size_t m = 0; m < m_wall.size(); m++) {
         m_wall[m]->initialize();
         if (m_wall[m]->kinetics(m_lr[m])) {
             nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies();
@@ -107,7 +107,7 @@ size_t Reactor::nSensParams()
         // determine the number of sensitivity parameters
         size_t m, ns;
         m_nsens = m_pnum.size();
-        for (m = 0; m < m_nwalls; m++) {
+        for (m = 0; m < m_wall.size(); m++) {
             ns = m_wall[m]->nSensParams(m_lr[m]);
             m_nsens_wall.push_back(ns);
             m_nsens += ns;
@@ -179,7 +179,7 @@ void Reactor::updateState(doublereal* y)
 void Reactor::updateSurfaceState(double* y)
 {
     size_t loc = 0;
-    for (size_t m = 0; m < m_nwalls; m++) {
+    for (size_t m = 0; m < m_wall.size(); m++) {
         SurfPhase* surf = m_wall[m]->surface(m_lr[m]);
         if (surf) {
             m_wall[m]->setCoverages(m_lr[m], y+loc);
@@ -230,31 +230,26 @@ void Reactor::evalEqs(doublereal time, doublereal* y,
         ydot[2] = 0.0;
     }
 
-    // add terms for open system
-    if (m_open) {
-        doublereal enthalpy = m_thermo->enthalpy_mass();
-
-        // outlets
-        for (size_t i = 0; i < m_nOutlets; i++) {
-            double mdot_out = m_outlet[i]->massFlowRate(time);
-            dmdt -= mdot_out; // mass flow out of system
-            if (m_energy) {
-                ydot[2] -= mdot_out * enthalpy;
-            }
+    // add terms for outlets
+    for (size_t i = 0; i < m_outlet.size(); i++) {
+        double mdot_out = m_outlet[i]->massFlowRate(time);
+        dmdt -= mdot_out; // mass flow out of system
+        if (m_energy) {
+            ydot[2] -= mdot_out * m_enthalpy;
         }
+    }
 
-        // inlets
-        for (size_t i = 0; i < m_nInlets; i++) {
-            double mdot_in = m_inlet[i]->massFlowRate(time);
-            dmdt += mdot_in; // mass flow into system
-            for (size_t n = 0; n < m_nsp; n++) {
-                double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
-                // flow of species into system and dilution by other species
-                dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
-            }
-            if (m_energy) {
-                ydot[2] += mdot_in * m_inlet[i]->enthalpy_mass();
-            }
+    // add terms for inlets
+    for (size_t i = 0; i < m_inlet.size(); i++) {
+        double mdot_in = m_inlet[i]->massFlowRate(time);
+        dmdt += mdot_in; // mass flow into system
+        for (size_t n = 0; n < m_nsp; n++) {
+            double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
+            // flow of species into system and dilution by other species
+            dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
+        }
+        if (m_energy) {
+            ydot[2] += mdot_in * m_inlet[i]->enthalpy_mass();
         }
     }
 
@@ -272,7 +267,7 @@ void Reactor::evalWalls(double t)
 {
     m_vdot = 0.0;
     m_Q = 0.0;
-    for (size_t i = 0; i < m_nwalls; i++) {
+    for (size_t i = 0; i < m_wall.size(); i++) {
         int lr = 1 - 2*m_lr[i];
         m_vdot += lr*m_wall[i]->vdot(t);
         m_Q += lr*m_wall[i]->Q(t);
@@ -287,7 +282,7 @@ double Reactor::evalSurfaces(double t, double* ydot)
     size_t loc = 0; // offset into ydot
     double mdot_surf = 0.0; // net mass flux from surface
 
-    for (size_t i = 0; i < m_nwalls; i++) {
+    for (size_t i = 0; i < m_wall.size(); i++) {
         Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
         SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
         if (surf && kin) {
@@ -332,7 +327,7 @@ std::vector<std::pair<void*, int> > Reactor::getSensitivityOrder() const
 {
     std::vector<std::pair<void*, int> > order;
     order.push_back(std::make_pair(const_cast<Reactor*>(this), 0));
-    for (size_t n = 0; n < m_nwalls; n++) {
+    for (size_t n = 0; n < m_wall.size(); n++) {
         if (m_nsens_wall[n]) {
             order.push_back(std::make_pair(m_wall[n], m_lr[n]));
         }
@@ -351,7 +346,7 @@ size_t Reactor::speciesIndex(const string& nm) const
     // check for a wall species
     size_t walloffset = 0, kp = 0;
     thermo_t* th;
-    for (size_t m = 0; m < m_nwalls; m++) {
+    for (size_t m = 0; m < m_wall.size(); m++) {
         if (m_wall[m]->kinetics(m_lr[m])) {
             kp = m_wall[m]->kinetics(m_lr[m])->reactionPhaseIndex();
             th = &m_wall[m]->kinetics(m_lr[m])->thermo(kp);
@@ -393,7 +388,7 @@ void Reactor::applySensitivity(double* params)
         m_kin->setMultiplier(m_pnum[n], mult*params[n]);
     }
     size_t ploc = npar;
-    for (size_t m = 0; m < m_nwalls; m++) {
+    for (size_t m = 0; m < m_wall.size(); m++) {
         if (m_nsens_wall[m] > 0) {
             m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc);
             ploc += m_nsens_wall[m];
@@ -413,7 +408,7 @@ void Reactor::resetSensitivity(double* params)
         m_kin->setMultiplier(m_pnum[n], mult/params[n]);
     }
     size_t ploc = npar;
-    for (size_t m = 0; m < m_nwalls; m++) {
+    for (size_t m = 0; m < m_wall.size(); m++) {
         if (m_nsens_wall[m] > 0) {
             m_wall[m]->resetSensitivityParameters(m_lr[m]);
             ploc += m_nsens_wall[m];
diff --git a/src/zeroD/ReactorBase.cpp b/src/zeroD/ReactorBase.cpp
index 9c60e43c7..5d2cef07f 100644
--- a/src/zeroD/ReactorBase.cpp
+++ b/src/zeroD/ReactorBase.cpp
@@ -6,7 +6,6 @@
 
 #include "cantera/zeroD/ReactorBase.h"
 #include "cantera/zeroD/FlowDevice.h"
-#include "cantera/zeroD/Wall.h"
 #include "cantera/zeroD/ReactorNet.h"
 
 using namespace std;
@@ -16,15 +15,10 @@ namespace Cantera
 ReactorBase::ReactorBase(const string& name) : m_nsp(0),
     m_thermo(0),
     m_vol(1.0),
-    m_vol0(1.0),
     m_init(false),
-    m_nInlets(0),
-    m_nOutlets(0),
-    m_open(false),
     m_enthalpy(0.0),
     m_intEnergy(0.0),
     m_pressure(0.0),
-    m_nwalls(0),
     m_net(0)
 {
     m_name = name;
@@ -54,15 +48,11 @@ void ReactorBase::syncState()
 void ReactorBase::addInlet(FlowDevice& inlet)
 {
     m_inlet.push_back(&inlet);
-    m_open = true;
-    m_nInlets++;
 }
 
 void ReactorBase::addOutlet(FlowDevice& outlet)
 {
     m_outlet.push_back(&outlet);
-    m_open = true;
-    m_nOutlets++;
 }
 
 void ReactorBase::addWall(Wall& w, int lr)
@@ -73,7 +63,6 @@ void ReactorBase::addWall(Wall& w, int lr)
     } else {
         m_lr.push_back(1);
     }
-    m_nwalls++;
 }
 
 Wall& ReactorBase::wall(size_t n)