diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index f8824a025..ab2121b9c 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -134,7 +134,7 @@ jobs:
             ~/nndc_hdf5
             ~/endf-b-vii.1
           key: ${{ runner.os }}-build-xs-cache
-      
+
       - name: before
         shell: bash
         run: $GITHUB_WORKSPACE/tools/ci/gha-before-script.sh
diff --git a/openmc/bounding_box.py b/openmc/bounding_box.py
index 4e1994283..daa6a318e 100644
--- a/openmc/bounding_box.py
+++ b/openmc/bounding_box.py
@@ -6,7 +6,7 @@ import numpy as np
 from .checkvalue import check_length
 
 
-class BoundingBox(tuple):
+class BoundingBox:
     """Axis-aligned bounding box.
 
     .. versionadded:: 0.13.4
@@ -35,17 +35,66 @@ class BoundingBox(tuple):
         The width of the x, y and z axis in [cm]
     """
 
-    def __new__(cls, lower_left: Iterable[float], upper_right: Iterable[float]):
+    def __init__(self, lower_left: Iterable[float], upper_right: Iterable[float]):
         check_length("lower_left", lower_left, 3, 3)
         check_length("upper_right", upper_right, 3, 3)
-        lower_left = np.array(lower_left, dtype=float)
-        upper_right = np.array(upper_right, dtype=float)
-        return tuple.__new__(cls, (lower_left, upper_right))
+        self._bounds = np.asarray([lower_left, upper_right], dtype=float)
 
     def __repr__(self) -> str:
         return "BoundingBox(lower_left={}, upper_right={})".format(
             tuple(self.lower_left), tuple(self.upper_right))
 
+    def __getitem__(self, key) -> np.ndarray:
+        return self._bounds[key]
+
+    def __len__(self):
+        return 2
+
+    def __setitem__(self, key, val):
+        self._bounds[key] = val
+
+    def __iand__(self, other: BoundingBox) -> BoundingBox:
+        """Updates the box be the intersection of itself and another box
+
+        Parameters
+        ----------
+        other : BoundingBox
+            The box used to resize this box
+
+        Returns
+        -------
+        An updated bounding box
+        """
+        self.lower_left = np.maximum(self.lower_left, other.lower_left)
+        self.upper_right = np.minimum(self.upper_right, other.upper_right)
+        return self
+
+    def __and__(self, other: BoundingBox) -> BoundingBox:
+        new = BoundingBox(*self)
+        new &= other
+        return new
+
+    def __ior__(self, other: BoundingBox) -> BoundingBox:
+        """Updates the box be the union of itself and another box
+
+        Parameters
+        ----------
+        other : BoundingBox
+            The box used to resize this box
+
+        Returns
+        -------
+        An updated bounding box
+        """
+        self.lower_left = np.minimum(self.lower_left, other.lower_left)
+        self.upper_right = np.maximum(self.upper_right, other.upper_right)
+        return self
+
+    def __or__(self, other: BoundingBox) -> BoundingBox:
+        new = BoundingBox(*self)
+        new |= other
+        return new
+
     @property
     def center(self) -> np.ndarray:
         return (self[0] + self[1]) / 2
@@ -54,10 +103,20 @@ class BoundingBox(tuple):
     def lower_left(self) -> np.ndarray:
         return self[0]
 
+    @lower_left.setter
+    def lower_left(self, llc):
+        check_length('lower_left', llc, 3, 3)
+        self[0] = llc
+
     @property
     def upper_right(self) -> np.ndarray:
         return self[1]
 
+    @upper_right.setter
+    def upper_right(self, urc):
+        check_length('upper_right', urc, 3, 3)
+        self[1] = urc
+
     @property
     def volume(self) -> float:
         return np.abs(np.prod(self[1] - self[0]))
@@ -89,29 +148,36 @@ class BoundingBox(tuple):
     def width(self):
         return self.upper_right - self.lower_left
 
-    def extend(self, padding_distance: float) -> BoundingBox:
-        """Returns an extended bounding box
+    def expand(self, padding_distance: float, inplace: bool = False) -> BoundingBox:
+        """Returns an expanded bounding box
 
         Parameters
         ----------
         padding_distance : float
             The distance to enlarge the bounding box by
+        inplace : bool
+            Whether or not to return a new BoundingBox instance or to modify the
+            current BoundingBox object.
 
         Returns
         -------
-        An enlarged bounding box
-
+        An expanded bounding box
         """
-        return BoundingBox(np.array(
-            [
-                self[0][0] - padding_distance,
-                self[0][1] - padding_distance,
-                self[0][2] - padding_distance
-            ]
-        ), np.array(
-            [
-                self[1][0] + padding_distance,
-                self[1][1] + padding_distance,
-                self[1][2] + padding_distance
-            ]
-        ))
+        if inplace:
+            self[0] -= padding_distance
+            self[1] += padding_distance
+            return self
+        else:
+            return BoundingBox(self[0] - padding_distance, self[1] + padding_distance)
+
+    @classmethod
+    def infinite(cls) -> BoundingBox:
+        """Create an infinite box. Useful as a starting point for determining
+           geometry bounds.
+
+        Returns
+        -------
+        An infinitely large bounding box.
+        """
+        infs = np.full((3,), np.inf)
+        return cls(-infs, infs)
diff --git a/openmc/region.py b/openmc/region.py
index 1d54bff8d..e0c05d674 100644
--- a/openmc/region.py
+++ b/openmc/region.py
@@ -414,13 +414,10 @@ class Intersection(Region, MutableSequence):
 
     @property
     def bounding_box(self):
-        lower_left = np.array([-np.inf, -np.inf, -np.inf])
-        upper_right = np.array([np.inf, np.inf, np.inf])
+        box = BoundingBox.infinite()
         for n in self:
-            lower_left_n, upper_right_n = n.bounding_box
-            lower_left[:] = np.maximum(lower_left, lower_left_n)
-            upper_right[:] = np.minimum(upper_right, upper_right_n)
-        return BoundingBox(lower_left, upper_right)
+            box &= n.bounding_box
+        return box
 
 
 class Union(Region, MutableSequence):
@@ -505,13 +502,11 @@ class Union(Region, MutableSequence):
 
     @property
     def bounding_box(self):
-        lower_left = np.array([np.inf, np.inf, np.inf])
-        upper_right = np.array([-np.inf, -np.inf, -np.inf])
+        bbox = BoundingBox(np.array([np.inf]*3),
+                           np.array([-np.inf]*3))
         for n in self:
-            lower_left_n, upper_right_n = n.bounding_box
-            lower_left[:] = np.minimum(lower_left, lower_left_n)
-            upper_right[:] = np.maximum(upper_right, upper_right_n)
-        return BoundingBox(lower_left, upper_right)
+            bbox |= n.bounding_box
+        return bbox
 
 
 class Complement(Region):
diff --git a/openmc/surface.py b/openmc/surface.py
index 6f19de5ca..dc1068b22 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -233,10 +233,7 @@ class Surface(IDManagerMixin, ABC):
             desired half-space
 
         """
-        return BoundingBox(
-            np.array([-np.inf, -np.inf, -np.inf]),
-            np.array([np.inf, np.inf, np.inf])
-        )
+        return BoundingBox.infinite()
 
     def clone(self, memo=None):
         """Create a copy of this surface with a new unique ID.
@@ -1204,13 +1201,9 @@ class Cylinder(QuadricMixin, Surface):
                   else -np.inf for xi, dxi in zip(self._origin, self._axis)]
             ur = [xi + r if np.isclose(dxi, 0., rtol=0., atol=self._atol)
                   else np.inf for xi, dxi in zip(self._origin, self._axis)]
-            return (np.array(ll), np.array(ur))
-
+            return BoundingBox(np.array(ll), np.array(ur))
         elif side == '+':
-            return BoundingBox(
-                np.array([-np.inf, -np.inf, -np.inf]),
-                np.array([np.inf, np.inf, np.inf])
-            )
+            return BoundingBox.infinite()
 
     def _get_base_coeffs(self):
         # Get x, y, z coordinates of two points
@@ -1377,10 +1370,7 @@ class XCylinder(QuadricMixin, Surface):
                 np.array([np.inf, self.y0 + self.r, self.z0 + self.r])
             )
         elif side == '+':
-            return BoundingBox(
-                np.array([-np.inf, -np.inf, -np.inf]),
-                np.array([np.inf, np.inf, np.inf])
-            )
+            return BoundingBox.infinite()
 
     def evaluate(self, point):
         y = point[1] - self.y0
@@ -1472,10 +1462,7 @@ class YCylinder(QuadricMixin, Surface):
                 np.array([self.x0 + self.r, np.inf, self.z0 + self.r])
             )
         elif side == '+':
-            return BoundingBox(
-                np.array([-np.inf, -np.inf, -np.inf]),
-                np.array([np.inf, np.inf, np.inf])
-            )
+            return BoundingBox.infinite()
 
     def evaluate(self, point):
         x = point[0] - self.x0
@@ -1567,10 +1554,7 @@ class ZCylinder(QuadricMixin, Surface):
                 np.array([self.x0 + self.r, self.y0 + self.r, np.inf])
             )
         elif side == '+':
-            return BoundingBox(
-                np.array([-np.inf, -np.inf, -np.inf]),
-                np.array([np.inf, np.inf, np.inf])
-            )
+            return BoundingBox.infinite()
 
     def evaluate(self, point):
         x = point[0] - self.x0
@@ -1661,10 +1645,7 @@ class Sphere(QuadricMixin, Surface):
                 np.array([self.x0 + self.r, self.y0 + self.r, self.z0 + self.r])
             )
         elif side == '+':
-            return BoundingBox(
-                np.array([-np.inf, -np.inf, -np.inf]),
-                np.array([np.inf, np.inf, np.inf])
-            )
+            return BoundingBox.infinite()
 
     def evaluate(self, point):
         x = point[0] - self.x0
@@ -2285,10 +2266,8 @@ class XTorus(TorusMixin, Surface):
                 np.array([x0 + b, y0 + a + c, z0 + a + c])
             )
         elif side == '+':
-            return BoundingBox(
-                np.array([-np.inf, -np.inf, -np.inf]),
-                np.array([np.inf, np.inf, np.inf])
-            )
+            return BoundingBox.infinite()
+
 
 class YTorus(TorusMixin, Surface):
     r"""A torus of the form :math:`(y - y_0)^2/B^2 + (\sqrt{(x - x_0)^2 + (z -
@@ -2360,10 +2339,7 @@ class YTorus(TorusMixin, Surface):
                 np.array([x0 + a + c, y0 + b, z0 + a + c])
             )
         elif side == '+':
-            return BoundingBox(
-                np.array([-np.inf, -np.inf, -np.inf]),
-                np.array([np.inf, np.inf, np.inf])
-            )
+            return BoundingBox.infinite()
 
 
 class ZTorus(TorusMixin, Surface):
@@ -2436,10 +2412,7 @@ class ZTorus(TorusMixin, Surface):
                 np.array([x0 + a + c, y0 + a + c, z0 + b])
             )
         elif side == '+':
-            return BoundingBox(
-                np.array([-np.inf, -np.inf, -np.inf]),
-                np.array([np.inf, np.inf, np.inf])
-            )
+            return BoundingBox.infinite()
 
 
 class Halfspace(Region):
diff --git a/openmc/universe.py b/openmc/universe.py
index cfa6634f8..44f95c1d2 100644
--- a/openmc/universe.py
+++ b/openmc/universe.py
@@ -234,8 +234,7 @@ class Universe(UniverseBase):
         if regions:
             return openmc.Union(regions).bounding_box
         else:
-            # Infinite bounding box
-            return openmc.Intersection([]).bounding_box
+            return openmc.BoundingBox.infinite()
 
     @classmethod
     def from_hdf5(cls, group, cells):
@@ -993,7 +992,7 @@ class DAGMCUniverse(UniverseBase):
         check_type('bounded type', bounded_type, str)
         check_value('bounded type', bounded_type, ('box', 'sphere'))
 
-        bbox = self.bounding_box.extend(padding_distance)
+        bbox = self.bounding_box.expand(padding_distance, True)
 
         if bounded_type == 'sphere':
             radius = np.linalg.norm(bbox.upper_right - bbox.center)
diff --git a/tests/unit_tests/test_bounding_box.py b/tests/unit_tests/test_bounding_box.py
index 1adb1976d..89e504271 100644
--- a/tests/unit_tests/test_bounding_box.py
+++ b/tests/unit_tests/test_bounding_box.py
@@ -81,3 +81,78 @@ def test_bounding_box_extents():
     assert test_bb_1.extent['xy'] == (-10., 1., -20., 2.)
     assert test_bb_1.extent['xz'] == (-10., 1., -30., 3.)
     assert test_bb_1.extent['yz'] == (-20., 2., -30., 3.)
+
+
+def test_bounding_box_methods():
+    test_bb = openmc.BoundingBox.infinite()
+
+    # check assignment operator
+    test_bb[0] = [-10, -11, -12]
+    test_bb[1] = [13, 14, 15]
+
+    assert all(test_bb[0] == [-10, -11, -12])
+    assert all(test_bb[1] == [13, 14, 15])
+
+    # check length and iteration
+    assert len(test_bb) == 2
+    ll, ur = test_bb
+    assert all(ll == [-10, -11, -12])
+    assert all(ur == [13, 14, 15])
+
+    # test expand/reduce methods
+    other_bb = openmc.BoundingBox([-5, -5, -50], [5, 50, 5])
+
+    reduced_bb = test_bb & other_bb
+
+    # inplace was False by default. BoundingBox.reduce should return a new object
+    assert test_bb is not reduced_bb
+
+    # the original bounding box should be unchanged
+    assert all(test_bb[0] == [-10, -11, -12])
+    assert all(test_bb[1] == [13, 14, 15])
+
+    assert all(reduced_bb[0] == [-5, -5, -12])
+    assert all(reduced_bb[1] == [5, 14, 5])
+
+    test_bb &= other_bb
+
+    assert all(test_bb[0] == [-5, -5, -12])
+    assert all(test_bb[1] == [5, 14, 5])
+
+    other_bb = openmc.BoundingBox([-50, -50, -1], [50, 1, 50])
+
+    expanded_bb = test_bb | other_bb
+
+    # inplace was False by default. BoundingBox.expand should return a new object
+    assert test_bb is not expanded_bb
+
+    # the original bounding box should be unchanged
+    assert all(test_bb[0] == [-5, -5, -12])
+    assert all(test_bb[1] == [5, 14, 5])
+
+    assert all(expanded_bb[0] == [-50, -50, -12])
+    assert all(expanded_bb[1] == [50, 14, 50])
+
+    test_bb |= other_bb
+
+    assert all(test_bb[0] == [-50, -50, -12])
+    assert all(test_bb[1] == [50, 14, 50])
+
+    extended_bbox = test_bb.expand(0.1)
+
+    assert extended_bbox is not test_bb
+
+    # the original bounding box should not be changed with inplace as False
+    assert all(test_bb[0] == [-50, -50, -12])
+    assert all(test_bb[1] == [50, 14, 50])
+
+    assert all(extended_bbox[0] == [-50.1, -50.1, -12.1])
+    assert all(extended_bbox[1] == [50.1, 14.1, 50.1])
+
+    extended_bbox = test_bb.expand(0.1, True)
+
+    # inplace was set to True. BoundingBox.reduce should return the same object
+    assert extended_bbox is test_bb
+
+    assert all(test_bb[0] == [-50.1, -50.1, -12.1])
+    assert all(test_bb[1] == [50.1, 14.1, 50.1])