diff --git a/src/base/SolutionArray.cpp b/src/base/SolutionArray.cpp
index 000c2328d1..4819aaa756 100644
--- a/src/base/SolutionArray.cpp
+++ b/src/base/SolutionArray.cpp
@@ -71,9 +71,19 @@ SolutionArray::SolutionArray(const SolutionArray& other,
     , m_extra(other.m_extra)
     , m_order(other.m_order)
     , m_shared(true)
-    , m_active(selected)
 {
     m_sol->thermo()->addSpeciesLock();
+    if (!other.m_shared) {
+        // direct slicing is possible
+        m_active = selected;
+    } else {
+        // slicing a previously sliced SolutionArray
+        m_active.clear();
+        m_active.reserve(selected.size());
+        for (auto loc : selected) {
+            m_active.push_back(other.m_active.at(loc));
+        }
+    }
     for (auto loc : m_active) {
         if (loc < 0 || loc >= (int)m_dataSize) {
             IndexError("SolutionArray::SolutionArray", "indices", loc, m_dataSize);
diff --git a/test/python/test_composite.py b/test/python/test_composite.py
index 260d09c2eb..58a40a0a6e 100644
--- a/test/python/test_composite.py
+++ b/test/python/test_composite.py
@@ -233,6 +233,24 @@ def test_from_state_array(self):
         arr = ct.SolutionArray(self.gas, states=states)
         assert arr.shape == (3, 5) # shape is based on numpy conversion
 
+    def test_slice_twice(self):
+        T_list = np.linspace(300, 1000, 8)
+        self.gas.TPX = T_list[0], ct.one_atm, {"H2": 1.}
+        arr = ct.SolutionArray(self.gas)
+        for T in T_list[1:]:
+            self.gas.TPX = T, ct.one_atm, {"H2": 1.}
+            arr.append(self.gas.state)
+        ix = 4
+        arr_trunc = arr[ix:]
+        assert arr_trunc.T[0] == arr.T[ix]
+        assert arr_trunc[0].T == arr.T[ix]
+        assert arr_trunc.T[-1] == arr.T[-1]
+        assert arr_trunc[-1].T == arr.T[-1]
+        assert (arr_trunc.T[:2] == arr.T[ix:ix+2]).all()
+        assert (arr_trunc[:2].T == arr.T[ix:ix+2]).all()
+        with self.assertRaises(IndexError):
+            arr_trunc[10]
+
     def test_from_state_numpy(self):
         states = np.array([[list(self.gas.state)] * 5] * 3)
         arr = ct.SolutionArray(self.gas, states=states)