diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst
index 97704e76d..a6c89be7c 100644
--- a/docs/source/pythonapi/model.rst
+++ b/docs/source/pythonapi/model.rst
@@ -25,7 +25,6 @@ Composite Surfaces
    :nosignatures:
    :template: myclass.rst
 
-   openmc.model.Box
    openmc.model.RectangularParallelepiped
    openmc.model.RightCircularCylinder
    openmc.model.XConeOneSided
diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 8c4b99ae4..7f444d4b1 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -156,85 +156,6 @@ class RectangularParallelepiped(CompositeSurface):
         return -self.xmin | +self.ymax | -self.ymin | +self.ymax | -self.zmin | +self.zmax
 
 
-class Box(CompositeSurface):
-    """Arbitrarily oriented orthogonal box
-
-    An arbitrarily oriented orthogonal box is composed of six planar surfaces.
-    This class acts as a proper surface, meaning that unary `+` and `-`
-    operators applied to it will produce a half-space. The negative side is
-    defined to be the region inside of the box.
-
-    .. versionadded:: 0.12
-
-    Parameters
-    ----------
-    v : iterable of float
-        Cartesian coordinates of a corner of the box
-    a1 : iterable of float
-        Vector of first side from the specified corner coordinates
-    a2 : iterable of float
-        Vector of second side from the specified corner coordinates
-    a3 : iterable of float
-        Vector of third side from the specified corner coordinates
-    **kwargs
-        Keyword arguments passed to underlying plane classes
-
-    Attributes
-    ----------
-    xmin, xmax : openmc.XPlane
-        Sides of the box
-    ymin, ymax : openmc.YPlane
-        Sides of the box
-    zmin, zmax : openmc.ZPlane
-        Sides of the box
-
-    """
-    _surface_names = ('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax')
-
-    def __init__(self, v, a1, a2, a3, **kwargs):
-        vx, vy, vz = v
-        a1x, a1y, a1z = a1
-        a2x, a2y, a2z = a2
-        a3x, a3y, a3z = a3
-
-        # Only support boxes with axis-aligned vectors
-        if any(x != 0.0 for x in (a1y, a1z, a2x, a2z, a3x, a3y)):
-            raise NotImplementedError('Box macrobody with non-axis-aligned '
-                                      'vector not supported.')
-
-        # Determine each side of the box
-        if a1x > 0:
-            xmin, xmax = vx, vx + a1x
-        else:
-            xmin, xmax = vx + a1x, vx
-        if a2y > 0:
-            ymin, ymax = vy, vy + a2y
-        else:
-            ymin, ymax = vy + a2y, vy
-        if a3z > 0:
-            zmin, zmax = vz, vz + a3z
-        else:
-            zmin, zmax = vz + a3z, vz
-
-        # Create surfaces
-        self.xmin = openmc.XPlane(x0=xmin, **kwargs)
-        self.xmax = openmc.XPlane(x0=xmax, **kwargs)
-        self.ymin = openmc.YPlane(y0=ymin, **kwargs)
-        self.ymax = openmc.YPlane(y0=ymax, **kwargs)
-        self.zmin = openmc.ZPlane(z0=zmin, **kwargs)
-        self.zmax = openmc.ZPlane(z0=zmax, **kwargs)
-
-    def __neg__(self):
-        return (+self.xmin & -self.xmax &
-                +self.ymin & -self.ymax &
-                +self.zmin & -self.zmax)
-
-    def __pos__(self):
-        return (-self.xmin | +self.xmax |
-                -self.ymin | +self.ymax |
-                -self.zmin | +self.zmax)
-
-
 class XConeOneSided(CompositeSurface):
     """One-sided cone parallel the x-axis
 
diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py
index 28ded03cb..a1a218311 100644
--- a/tests/unit_tests/test_surface_composite.py
+++ b/tests/unit_tests/test_surface_composite.py
@@ -5,63 +5,6 @@ import openmc
 import pytest
 
 
-def test_box():
-    x, y, z = 1.0, -2.5, 3.0
-    dx, dy, dz = 4.0, 3.0, 6.0
-    s = openmc.model.Box((x, y, z), (dx, 0, 0), (0, dy, 0), (0, 0, dz))
-    assert isinstance(s.xmin, openmc.XPlane)
-    assert isinstance(s.xmax, openmc.XPlane)
-    assert isinstance(s.ymin, openmc.YPlane)
-    assert isinstance(s.ymax, openmc.YPlane)
-    assert isinstance(s.zmin, openmc.ZPlane)
-    assert isinstance(s.zmax, openmc.ZPlane)
-
-    # Make sure boundary condition propagates
-    s.boundary_type = 'reflective'
-    assert s.boundary_type == 'reflective'
-    for axis in 'xyz':
-        assert getattr(s, '{}min'.format(axis)).boundary_type == 'reflective'
-        assert getattr(s, '{}max'.format(axis)).boundary_type == 'reflective'
-
-    # Check bounding box
-    ll, ur = (+s).bounding_box
-    assert np.all(np.isinf(ll))
-    assert np.all(np.isinf(ur))
-    ll, ur = (-s).bounding_box
-    assert ur == pytest.approx((x + dx, y + dy, z + dz))
-    assert ll == pytest.approx((x, y, z))
-
-    # __contains__ on associated half-spaces
-    assert (x - 1, 0, 0) in +s
-    assert (x - 1, 0, 0) not in -s
-    assert (x + dx/2, y + dy/2, z + dz/2) in -s
-    assert (x + dx/2, y + dy/2, z + dz/2) not in +s
-
-    # evaluate method
-    with pytest.raises(NotImplementedError):
-        s.evaluate((0., 0., 0.))
-
-    # translate method
-    t = uniform(-5.0, 5.0)
-    s_t = s.translate((t, t, t))
-    ll_t, ur_t = (-s_t).bounding_box
-    assert ur_t == pytest.approx(ur + t)
-    assert ll_t == pytest.approx(ll + t)
-
-    # rotate method
-    s_rot = s.rotate((0., 0., 90.))
-    ll, ur = (-s_rot).bounding_box
-    assert ll == pytest.approx((-(y + dy), x, z))
-    assert ur == pytest.approx((-y, x + dx, z + dz))
-
-    # Make sure repr works
-    repr(s)
-
-    # non-axis-aligned vectors not implemented
-    with pytest.raises(NotImplementedError):
-        openmc.model.Box((0, 0, 0), (1, 1, 1), (1, 1, 1), (1, 1, 1))
-
-
 def test_rectangular_parallelepiped():
     xmin = uniform(-5., 5.)
     xmax = xmin + uniform(0., 5.)