diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst
index 72d2534ea..d379fb4bb 100644
--- a/docs/source/pythonapi/base.rst
+++ b/docs/source/pythonapi/base.rst
@@ -87,6 +87,7 @@ Building geometry
    openmc.Intersection
    openmc.Union
    openmc.Complement
+   openmc.BoundingBox
    openmc.Cell
    openmc.Universe
    openmc.DAGMCUniverse
diff --git a/openmc/__init__.py b/openmc/__init__.py
index ac9e02e16..5bb150fb0 100644
--- a/openmc/__init__.py
+++ b/openmc/__init__.py
@@ -1,4 +1,5 @@
 from openmc.arithmetic import *
+from openmc.bounding_box import *
 from openmc.cell import *
 from openmc.checkvalue import *
 from openmc.mesh import *
diff --git a/openmc/bounding_box.py b/openmc/bounding_box.py
new file mode 100644
index 000000000..d2a499092
--- /dev/null
+++ b/openmc/bounding_box.py
@@ -0,0 +1,57 @@
+from typing import Iterable
+
+import numpy as np
+
+from .checkvalue import check_length
+
+
+class BoundingBox(tuple):
+    """Axis-aligned bounding box.
+
+    .. versionadded:: 0.13.4
+
+    Parameters
+    ----------
+    lower_left : iterable of float
+        The x, y, z coordinates of the lower left corner of the bounding box in [cm]
+    upper_right : iterable of float
+        The x, y, z coordinates of the upper right corner of the bounding box in [cm]
+
+    Attributes
+    ----------
+    center : numpy.ndarray
+        x, y, z coordinates of the center of the bounding box in [cm]
+    lower_left : numpy.ndarray
+        The x, y, z coordinates of the lower left corner of the bounding box in [cm]
+    upper_right : numpy.ndarray
+        The x, y, z coordinates of the upper right corner of the bounding box in [cm]
+    volume : float
+        The volume of the bounding box in [cm^3]
+    """
+
+    def __new__(cls, 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))
+
+    def __repr__(self) -> str:
+        return "BoundingBox(lower_left={}, upper_right={})".format(
+            tuple(self.lower_left), tuple(self.upper_right))
+
+    @property
+    def center(self) -> np.ndarray:
+        return (self[0] + self[1]) / 2
+
+    @property
+    def lower_left(self) -> np.ndarray:
+        return self[0]
+
+    @property
+    def upper_right(self) -> np.ndarray:
+        return self[1]
+
+    @property
+    def volume(self) -> float:
+        return np.abs(np.prod(self[1] - self[0]))
diff --git a/openmc/cell.py b/openmc/cell.py
index a8e1178f4..67735ef5e 100644
--- a/openmc/cell.py
+++ b/openmc/cell.py
@@ -13,6 +13,7 @@ from ._xml import get_text
 from .mixin import IDManagerMixin
 from .region import Region, Complement
 from .surface import Halfspace
+from .bounding_box import BoundingBox
 
 
 class Cell(IDManagerMixin):
@@ -91,6 +92,8 @@ class Cell(IDManagerMixin):
         fill. For example, {'U235': 1.0e22, 'U238': 5.0e22, ...}.
 
         .. versionadded:: 0.12
+    bounding_box : openmc.BoundingBox
+        Axis-aligned bounding box of the cell
 
     """
 
@@ -249,8 +252,8 @@ class Cell(IDManagerMixin):
         if self.region is not None:
             return self.region.bounding_box
         else:
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([np.inf, np.inf, np.inf]))
+            return BoundingBox(np.array([-np.inf, -np.inf, -np.inf]),
+                               np.array([np.inf, np.inf, np.inf]))
 
     @property
     def num_instances(self):
diff --git a/openmc/geometry.py b/openmc/geometry.py
index ee4a58a3b..dbfce731f 100644
--- a/openmc/geometry.py
+++ b/openmc/geometry.py
@@ -28,7 +28,7 @@ class Geometry:
     ----------
     root_universe : openmc.UniverseBase
         Root universe which contains all others
-    bounding_box : 2-tuple of numpy.array
+    bounding_box : openmc.BoundingBox
         Lower-left and upper-right coordinates of an axis-aligned bounding box
         of the universe.
     merge_surfaces : bool
diff --git a/openmc/lib/cell.py b/openmc/lib/cell.py
index b53a6be7a..807fe73c7 100644
--- a/openmc/lib/cell.py
+++ b/openmc/lib/cell.py
@@ -11,6 +11,8 @@ from . import _dll
 from .core import _FortranObjectWithID
 from .error import _error_handler
 from .material import Material
+from ..bounding_box import BoundingBox
+
 
 __all__ = ['Cell', 'cells']
 
@@ -101,8 +103,8 @@ class Cell(_FortranObjectWithID):
         Name of the cell
     num_instances : int
         Number of unique cell instances
-    bounding_box : 2-tuple of numpy.ndarray
-        Lower-left and upper-right coordinates of bounding box
+    bounding_box : openmc.BoundingBox
+        Axis-aligned bounding box of the cell
     translation : Iterable of float
         3-D coordinates of the translation vector
     rotation : Iterable of float
@@ -289,7 +291,7 @@ class Cell(_FortranObjectWithID):
         llc[llc == -inf] = -np.inf
         urc[urc == -inf] = -np.inf
 
-        return llc, urc
+        return BoundingBox(llc, urc)
 
 
 class _CellMapping(Mapping):
diff --git a/openmc/region.py b/openmc/region.py
index 5d0e68020..af9df1ba1 100644
--- a/openmc/region.py
+++ b/openmc/region.py
@@ -6,6 +6,7 @@ from copy import deepcopy
 import numpy as np
 
 from .checkvalue import check_type
+from .bounding_box import BoundingBox
 
 
 class Region(ABC):
@@ -355,7 +356,7 @@ class Intersection(Region, MutableSequence):
 
     Attributes
     ----------
-    bounding_box : tuple of numpy.array
+    bounding_box : openmc.BoundingBox
         Lower-left and upper-right coordinates of an axis-aligned bounding box
 
     """
@@ -418,7 +419,7 @@ class Intersection(Region, MutableSequence):
             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 lower_left, upper_right
+        return BoundingBox(lower_left, upper_right)
 
 
 class Union(Region, MutableSequence):
@@ -443,7 +444,7 @@ class Union(Region, MutableSequence):
 
     Attributes
     ----------
-    bounding_box : 2-tuple of numpy.array
+    bounding_box : openmc.BoundingBox
         Lower-left and upper-right coordinates of an axis-aligned bounding box
 
     """
@@ -506,7 +507,7 @@ class Union(Region, MutableSequence):
             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 lower_left, upper_right
+        return BoundingBox(lower_left, upper_right)
 
 
 class Complement(Region):
@@ -533,7 +534,7 @@ class Complement(Region):
     ----------
     node : openmc.Region
         Regions to take the complement of
-    bounding_box : tuple of numpy.array
+    bounding_box : openmc.BoundingBox
         Lower-left and upper-right coordinates of an axis-aligned bounding box
 
     """
diff --git a/openmc/surface.py b/openmc/surface.py
index 984bc7a5c..5d13b39ff 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -12,6 +12,7 @@ import numpy as np
 from .checkvalue import check_type, check_value, check_length
 from .mixin import IDManagerMixin, IDWarning
 from .region import Region, Intersection, Union
+from .bounding_box import BoundingBox
 
 
 _BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic', 'white']
@@ -233,8 +234,10 @@ class Surface(IDManagerMixin, ABC):
             desired half-space
 
         """
-        return (np.array([-np.inf, -np.inf, -np.inf]),
-                np.array([np.inf, np.inf, np.inf]))
+        return BoundingBox(
+            np.array([-np.inf, -np.inf, -np.inf]),
+            np.array([np.inf, np.inf, np.inf])
+        )
 
     def clone(self, memo=None):
         """Create a copy of this surface with a new unique ID.
@@ -539,7 +542,7 @@ class PlaneMixin:
                 else:
                     ur = np.array([v if not np.isnan(v) else np.inf for v in vals])
 
-        return (ll, ur)
+        return BoundingBox(ll, ur)
 
     def evaluate(self, point):
         """Evaluate the surface equation at a given point.
@@ -1205,8 +1208,10 @@ class Cylinder(QuadricMixin, Surface):
             return (np.array(ll), np.array(ur))
 
         elif side == '+':
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([np.inf, np.inf, np.inf]))
+            return BoundingBox(
+                np.array([-np.inf, -np.inf, -np.inf]),
+                np.array([np.inf, np.inf, np.inf])
+            )
 
     def _get_base_coeffs(self):
         # Get x, y, z coordinates of two points
@@ -1368,11 +1373,15 @@ class XCylinder(QuadricMixin, Surface):
 
     def bounding_box(self, side):
         if side == '-':
-            return (np.array([-np.inf, self.y0 - self.r, self.z0 - self.r]),
-                    np.array([np.inf, self.y0 + self.r, self.z0 + self.r]))
+            return BoundingBox(
+                np.array([-np.inf, self.y0 - self.r, self.z0 - self.r]),
+                np.array([np.inf, self.y0 + self.r, self.z0 + self.r])
+            )
         elif side == '+':
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([np.inf, np.inf, np.inf]))
+            return BoundingBox(
+                np.array([-np.inf, -np.inf, -np.inf]),
+                np.array([np.inf, np.inf, np.inf])
+            )
 
     def evaluate(self, point):
         y = point[1] - self.y0
@@ -1459,11 +1468,15 @@ class YCylinder(QuadricMixin, Surface):
 
     def bounding_box(self, side):
         if side == '-':
-            return (np.array([self.x0 - self.r, -np.inf, self.z0 - self.r]),
-                    np.array([self.x0 + self.r, np.inf, self.z0 + self.r]))
+            return BoundingBox(
+                np.array([self.x0 - self.r, -np.inf, self.z0 - self.r]),
+                np.array([self.x0 + self.r, np.inf, self.z0 + self.r])
+            )
         elif side == '+':
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([np.inf, np.inf, np.inf]))
+            return BoundingBox(
+                np.array([-np.inf, -np.inf, -np.inf]),
+                np.array([np.inf, np.inf, np.inf])
+            )
 
     def evaluate(self, point):
         x = point[0] - self.x0
@@ -1550,11 +1563,15 @@ class ZCylinder(QuadricMixin, Surface):
 
     def bounding_box(self, side):
         if side == '-':
-            return (np.array([self.x0 - self.r, self.y0 - self.r, -np.inf]),
-                    np.array([self.x0 + self.r, self.y0 + self.r, np.inf]))
+            return BoundingBox(
+                np.array([self.x0 - self.r, self.y0 - self.r, -np.inf]),
+                np.array([self.x0 + self.r, self.y0 + self.r, np.inf])
+            )
         elif side == '+':
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([np.inf, np.inf, np.inf]))
+            return BoundingBox(
+                np.array([-np.inf, -np.inf, -np.inf]),
+                np.array([np.inf, np.inf, np.inf])
+            )
 
     def evaluate(self, point):
         x = point[0] - self.x0
@@ -1640,13 +1657,15 @@ class Sphere(QuadricMixin, Surface):
 
     def bounding_box(self, side):
         if side == '-':
-            return (np.array([self.x0 - self.r, self.y0 - self.r,
-                              self.z0 - self.r]),
-                    np.array([self.x0 + self.r, self.y0 + self.r,
-                              self.z0 + self.r]))
+            return BoundingBox(
+                np.array([self.x0 - self.r, self.y0 - self.r, self.z0 - self.r]),
+                np.array([self.x0 + self.r, self.y0 + self.r, self.z0 + self.r])
+            )
         elif side == '+':
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([np.inf, np.inf, np.inf]))
+            return BoundingBox(
+                np.array([-np.inf, -np.inf, -np.inf]),
+                np.array([np.inf, np.inf, np.inf])
+            )
 
     def evaluate(self, point):
         x = point[0] - self.x0
@@ -2262,11 +2281,15 @@ class XTorus(TorusMixin, Surface):
         x0, y0, z0 = self.x0, self.y0, self.z0
         a, b, c = self.a, self.b, self.c
         if side == '-':
-            return (np.array([x0 - b, y0 - a - c, z0 - a - c]),
-                    np.array([x0 + b, y0 + a + c, z0 + a + c]))
+            return BoundingBox(
+                np.array([x0 - b, y0 - a - c, z0 - a - c]),
+                np.array([x0 + b, y0 + a + c, z0 + a + c])
+            )
         elif side == '+':
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([np.inf, np.inf, np.inf]))
+            return BoundingBox(
+                np.array([-np.inf, -np.inf, -np.inf]),
+                np.array([np.inf, np.inf, np.inf])
+            )
 
 class YTorus(TorusMixin, Surface):
     r"""A torus of the form :math:`(y - y_0)^2/B^2 + (\sqrt{(x - x_0)^2 + (z -
@@ -2333,11 +2356,16 @@ class YTorus(TorusMixin, Surface):
         x0, y0, z0 = self.x0, self.y0, self.z0
         a, b, c = self.a, self.b, self.c
         if side == '-':
-            return (np.array([x0 - a - c, y0 - b, z0 - a - c]),
-                    np.array([x0 + a + c, y0 + b, z0 + a + c]))
+            return BoundingBox(
+                np.array([x0 - a - c, y0 - b, z0 - a - c]),
+                np.array([x0 + a + c, y0 + b, z0 + a + c])
+            )
         elif side == '+':
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([np.inf, np.inf, np.inf]))
+            return BoundingBox(
+                np.array([-np.inf, -np.inf, -np.inf]),
+                np.array([np.inf, np.inf, np.inf])
+            )
+
 
 class ZTorus(TorusMixin, Surface):
     r"""A torus of the form :math:`(z - z_0)^2/B^2 + (\sqrt{(x - x_0)^2 + (y -
@@ -2404,11 +2432,15 @@ class ZTorus(TorusMixin, Surface):
         x0, y0, z0 = self.x0, self.y0, self.z0
         a, b, c = self.a, self.b, self.c
         if side == '-':
-            return (np.array([x0 - a - c, y0 - a - c, z0 - b]),
-                    np.array([x0 + a + c, y0 + a + c, z0 + b]))
+            return BoundingBox(
+                np.array([x0 - a - c, y0 - a - c, z0 - b]),
+                np.array([x0 + a + c, y0 + a + c, z0 + b])
+            )
         elif side == '+':
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([np.inf, np.inf, np.inf]))
+            return BoundingBox(
+                np.array([-np.inf, -np.inf, -np.inf]),
+                np.array([np.inf, np.inf, np.inf])
+            )
 
 
 class Halfspace(Region):
@@ -2443,7 +2475,7 @@ class Halfspace(Region):
         Surface which divides Euclidean space.
     side : {'+', '-'}
         Indicates whether the positive or negative half-space is used.
-    bounding_box : tuple of numpy.ndarray
+    bounding_box : openmc.BoundingBox
         Lower-left and upper-right coordinates of an axis-aligned bounding box
 
     """
diff --git a/openmc/universe.py b/openmc/universe.py
index 67296dda8..ab9858f68 100644
--- a/openmc/universe.py
+++ b/openmc/universe.py
@@ -205,7 +205,7 @@ class Universe(UniverseBase):
         Volume of the universe in cm^3. This can either be set manually or
         calculated in a stochastic volume calculation and added via the
         :meth:`Universe.add_volume_information` method.
-    bounding_box : 2-tuple of numpy.array
+    bounding_box : openmc.BoundingBox
         Lower-left and upper-right coordinates of an axis-aligned bounding box
         of the universe.
 
@@ -733,7 +733,7 @@ class DAGMCUniverse(UniverseBase):
     auto_mat_ids : bool
         Set IDs automatically on initialization (True)  or report overlaps in ID
         space between OpenMC and UWUW materials (False)
-    bounding_box : 2-tuple of numpy.array
+    bounding_box : openmc.BoundingBox
         Lower-left and upper-right coordinates of an axis-aligned bounding box
         of the universe.
 
diff --git a/tests/unit_tests/test_bounding_box.py b/tests/unit_tests/test_bounding_box.py
new file mode 100644
index 000000000..f053bb2a1
--- /dev/null
+++ b/tests/unit_tests/test_bounding_box.py
@@ -0,0 +1,77 @@
+import numpy as np
+import openmc
+import pytest
+
+
+test_bb_1 = openmc.BoundingBox((-10.0, -20.0, -30.0), (1.0, 2.0, 3.0))
+test_bb_2 = openmc.BoundingBox((1.0, 2.0, 3.0), (11.0, 22.0, 33.0))
+test_bb_3 = openmc.BoundingBox((-10.0, -20.0, -30.0), (-1.0, -2.0, -3.0))
+
+
+@pytest.mark.parametrize(
+    "bb, expected",
+    [
+        (test_bb_1, 7986),  # 11 * 22 * 33
+        (test_bb_2, 6000),  # 10 * 20 * 30
+        (test_bb_3, 4374),  # 9 * 18 * 27
+    ],
+)
+def test_bounding_box_volume(bb, expected):
+    assert bb.volume == expected
+
+
+@pytest.mark.parametrize(
+    "bb, expected",
+    [
+        (test_bb_1, np.array([-10.0, -20.0, -30.0])),
+        (test_bb_2, np.array([1.0, 2.0, 3.0])),
+        (test_bb_3, np.array([-10.0, -20.0, -30.0])),
+    ],
+)
+def test_bounding_lower_left(bb, expected):
+    assert np.array_equiv(expected, bb.lower_left)
+
+
+@pytest.mark.parametrize(
+    "bb, expected",
+    [
+        (test_bb_1, np.array([1.0, 2.0, 3.0])),
+        (test_bb_2, np.array([11.0, 22.0, 33.0])),
+        (test_bb_3, np.array([-1.0, -2.0, -3.0])),
+    ],
+)
+def test_bounding_upper_right(bb, expected):
+    assert np.array_equiv(expected, bb.upper_right)
+
+
+@pytest.mark.parametrize(
+    "bb, expected",
+    [
+        (test_bb_1, np.array([-4.5, -9.0, -13.5])),
+        (test_bb_2, np.array([6.0, 12.0, 18.0])),
+        (test_bb_3, np.array([-5.5, -11.0, -16.5])),
+    ],
+)
+def test_bounding_box_center(bb, expected):
+    assert np.array_equiv(expected, bb.center)
+
+
+def test_bounding_box_input_checking():
+    # checks that only passing lower_left is not accepted
+    with pytest.raises(TypeError):
+        openmc.BoundingBox((-10, -20, -3))
+    # checks that a tuple with three entry is not accepted
+    with pytest.raises(TypeError):
+        openmc.BoundingBox((-1, -2, -3), (-1, -2, -3), (-1, -2, -3))
+    # checks that a numpy array with two entries is not accepted
+    with pytest.raises(ValueError):
+        openmc.BoundingBox(np.array([-10, -30]), np.array([1, 2, 3]))
+    # checks that a numpy array with two entries is not accepted
+    with pytest.raises(ValueError):
+        openmc.BoundingBox(np.array([-10, -20, -30]), np.array([1, 3]))
+    # checks that a numpy array with four entries is not accepted
+    with pytest.raises(ValueError):
+        openmc.BoundingBox(np.array([-10, -20, -3, -4]), np.array([1, 2, 3]))
+    # checks that a numpy array with four entries is not accepted
+    with pytest.raises(ValueError):
+        openmc.BoundingBox(np.array([-10, -20, -4]), np.array([1, 2, 3, 4]))