diff --git a/.gitmodules b/.gitmodules
index 4dde5fd5ee5..f84d09bb1f9 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,9 +1,6 @@
 [submodule "vendor/pugixml"]
 	path = vendor/pugixml
 	url = https://github.com/zeux/pugixml.git
-[submodule "vendor/gsl-lite"]
-	path = vendor/gsl-lite
-	url = https://github.com/martinmoene/gsl-lite.git
 [submodule "vendor/xtensor"]
 	path = vendor/xtensor
 	url = https://github.com/xtensor-stack/xtensor.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 575e45373ae..9256204f3af 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -274,19 +274,6 @@ if (NOT xtensor_FOUND)
   add_subdirectory(vendor/xtensor)
 endif()
 
-#===============================================================================
-# GSL header-only library
-#===============================================================================
-
-find_package_write_status(gsl-lite)
-if (NOT gsl-lite_FOUND)
-  add_subdirectory(vendor/gsl-lite)
-
-  # Make sure contract violations throw exceptions
-  target_compile_definitions(gsl-lite-v1 INTERFACE GSL_THROW_ON_CONTRACT_VIOLATION)
-  target_compile_definitions(gsl-lite-v1 INTERFACE gsl_CONFIG_ALLOWS_NONSTRICT_SPAN_COMPARISON=1)
-endif()
-
 #===============================================================================
 # Catch2 library
 #===============================================================================
@@ -499,7 +486,7 @@ endif()
 # target_link_libraries treats any arguments starting with - but not -l as
 # linker flags. Thus, we can pass both linker flags and libraries together.
 target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}
-                      xtensor gsl::gsl-lite-v1 fmt::fmt ${CMAKE_DL_LIBS})
+                      xtensor fmt::fmt ${CMAKE_DL_LIBS})
 
 if(TARGET pugixml::pugixml)
   target_link_libraries(libopenmc pugixml::pugixml)
diff --git a/MANIFEST.in b/MANIFEST.in
index b73218af0dc..cdc7e2abcf0 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -43,6 +43,5 @@ recursive-include vendor *.hh
 recursive-include vendor *.hpp
 recursive-include vendor *.pc.in
 recursive-include vendor *.natvis
-include vendor/gsl-lite/include/gsl/gsl
 prune docs/build
 prune docs/source/pythonapi/generated/
diff --git a/cmake/OpenMCConfig.cmake.in b/cmake/OpenMCConfig.cmake.in
index b3b901de427..f90721b9319 100644
--- a/cmake/OpenMCConfig.cmake.in
+++ b/cmake/OpenMCConfig.cmake.in
@@ -1,7 +1,6 @@
 get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY)
 
 find_package(fmt REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../fmt)
-find_package(gsl-lite REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../gsl-lite)
 find_package(pugixml REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../pugixml)
 find_package(xtl REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtl)
 find_package(xtensor REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtensor)
@@ -41,4 +40,4 @@ endif()
 
 if(@OPENMC_USE_UWUW@ AND NOT ${DAGMC_BUILD_UWUW})
   message(FATAL_ERROR "UWUW is enabled in OpenMC but the DAGMC installation discovered was not configured with UWUW.")
-endif()
\ No newline at end of file
+endif()
diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 70843140bad..516a3236192 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -10,7 +10,6 @@
 
 #include "hdf5.h"
 #include "pugixml.hpp"
-#include <gsl/gsl-lite.hpp>
 
 #include "openmc/bounding_box.h"
 #include "openmc/constants.h"
@@ -129,7 +128,7 @@ class Region {
   void add_precedence();
 
   //! Add parenthesis to enforce precedence
-  gsl::index add_parentheses(gsl::index start);
+  int64_t add_parentheses(int64_t start);
 
   //! Remove complement operators from the expression
   void remove_complement_ops();
@@ -415,8 +414,8 @@ struct CellInstance {
     return index_cell == other.index_cell && instance == other.instance;
   }
 
-  gsl::index index_cell;
-  gsl::index instance;
+  int64_t index_cell;
+  int64_t instance;
 };
 
 //! Structure necessary for inserting CellInstance into hashed STL data
diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index e70c803de2f..854cf7d7719 100644
--- a/include/openmc/distribution.h
+++ b/include/openmc/distribution.h
@@ -7,10 +7,10 @@
 #include <cstddef> // for size_t
 
 #include "pugixml.hpp"
-#include <gsl/gsl-lite.hpp>
 
 #include "openmc/constants.h"
 #include "openmc/memory.h" // for unique_ptr
+#include "openmc/span.h"
 #include "openmc/vector.h" // for vector
 
 namespace openmc {
@@ -44,9 +44,9 @@ class DiscreteIndex {
 public:
   DiscreteIndex() {};
   DiscreteIndex(pugi::xml_node node);
-  DiscreteIndex(gsl::span<const double> p);
+  DiscreteIndex(span<const double> p);
 
-  void assign(gsl::span<const double> p);
+  void assign(span<const double> p);
 
   //! Sample a value from the distribution
   //! \param seed Pseudorandom number seed pointer
diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h
index 8ff766d1333..9c3bc743ffb 100644
--- a/include/openmc/distribution_spatial.h
+++ b/include/openmc/distribution_spatial.h
@@ -6,6 +6,7 @@
 #include "openmc/distribution.h"
 #include "openmc/mesh.h"
 #include "openmc/position.h"
+#include "openmc/span.h"
 
 namespace openmc {
 
@@ -104,7 +105,7 @@ class SphericalIndependent : public SpatialDistribution {
 class MeshSpatial : public SpatialDistribution {
 public:
   explicit MeshSpatial(pugi::xml_node node);
-  explicit MeshSpatial(int32_t mesh_id, gsl::span<const double> strengths);
+  explicit MeshSpatial(int32_t mesh_id, span<const double> strengths);
 
   //! Sample a position from the distribution
   //! \param seed Pseudorandom number seed pointer
@@ -144,7 +145,7 @@ class PointCloud : public SpatialDistribution {
 public:
   explicit PointCloud(pugi::xml_node node);
   explicit PointCloud(
-    std::vector<Position> point_cloud, gsl::span<const double> strengths);
+    std::vector<Position> point_cloud, span<const double> strengths);
 
   //! Sample a position from the distribution
   //! \param seed Pseudorandom number seed pointer
diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h
index db501f71b66..31ae8b0e3f7 100644
--- a/include/openmc/interpolate.h
+++ b/include/openmc/interpolate.h
@@ -4,10 +4,9 @@
 #include <cmath>
 #include <vector>
 
-#include <gsl/gsl-lite.hpp>
-
 #include "openmc/error.h"
 #include "openmc/search.h"
+#include "openmc/span.h"
 
 namespace openmc {
 
@@ -36,8 +35,8 @@ inline double interpolate_log_log(
   return y0 * std::exp(f * std::log(y1 / y0));
 }
 
-inline double interpolate_lagrangian(gsl::span<const double> xs,
-  gsl::span<const double> ys, int idx, double x, int order)
+inline double interpolate_lagrangian(
+  span<const double> xs, span<const double> ys, int idx, double x, int order)
 {
   double output {0.0};
 
@@ -56,9 +55,8 @@ inline double interpolate_lagrangian(gsl::span<const double> xs,
   return output;
 }
 
-inline double interpolate(gsl::span<const double> xs,
-  gsl::span<const double> ys, double x,
-  Interpolation i = Interpolation::lin_lin)
+inline double interpolate(span<const double> xs, span<const double> ys,
+  double x, Interpolation i = Interpolation::lin_lin)
 {
   int idx = lower_bound_index(xs.begin(), xs.end(), x);
 
diff --git a/include/openmc/material.h b/include/openmc/material.h
index 9235be356f6..31fab2ae273 100644
--- a/include/openmc/material.h
+++ b/include/openmc/material.h
@@ -4,9 +4,9 @@
 #include <string>
 #include <unordered_map>
 
+#include "openmc/span.h"
 #include "pugixml.hpp"
 #include "xtensor/xtensor.hpp"
-#include <gsl/gsl-lite.hpp>
 #include <hdf5.h>
 
 #include "openmc/bremsstrahlung.h"
@@ -118,21 +118,21 @@ class Material {
   //
   //! \param[in] density Density value
   //! \param[in] units Units of density
-  void set_density(double density, gsl::cstring_span units);
+  void set_density(double density, const std::string& units);
 
   //! Set temperature of the material
   void set_temperature(double temperature) { temperature_ = temperature; };
 
   //! Get nuclides in material
   //! \return Indices into the global nuclides vector
-  gsl::span<const int> nuclides() const
+  span<const int> nuclides() const
   {
     return {nuclide_.data(), nuclide_.size()};
   }
 
   //! Get densities of each nuclide in material
   //! \return Densities in [atom/b-cm]
-  gsl::span<const double> densities() const
+  span<const double> densities() const
   {
     return {atom_density_.data(), atom_density_.size()};
   }
@@ -210,7 +210,7 @@ class Material {
 
   //----------------------------------------------------------------------------
   // Private data members
-  gsl::index index_;
+  int64_t index_;
 
   bool depletable_ {false}; //!< Is the material depletable?
   bool fissionable_ {
diff --git a/include/openmc/mcpl_interface.h b/include/openmc/mcpl_interface.h
index 1f0c94d6dec..780f7e86321 100644
--- a/include/openmc/mcpl_interface.h
+++ b/include/openmc/mcpl_interface.h
@@ -2,10 +2,9 @@
 #define OPENMC_MCPL_INTERFACE_H
 
 #include "openmc/particle_data.h"
+#include "openmc/span.h"
 #include "openmc/vector.h"
 
-#include <gsl/gsl-lite.hpp>
-
 #include <string>
 
 namespace openmc {
@@ -36,7 +35,7 @@ vector<SourceSite> mcpl_source_sites(std::string path);
 //!                         calculate_parallel_index_vector on
 //!                         source_bank.size().
 void write_mcpl_source_point(const char* filename,
-  gsl::span<SourceSite> source_bank, vector<int64_t> const& bank_index);
+  span<const SourceSite> source_bank, const vector<int64_t>& bank_index);
 } // namespace openmc
 
 #endif // OPENMC_MCPL_INTERFACE_H
diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 5edb03003f9..1bd5160f193 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -9,13 +9,13 @@
 #include "hdf5.h"
 #include "pugixml.hpp"
 #include "xtensor/xtensor.hpp"
-#include <gsl/gsl-lite.hpp>
 
 #include "openmc/bounding_box.h"
 #include "openmc/error.h"
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/particle.h"
 #include "openmc/position.h"
+#include "openmc/span.h"
 #include "openmc/vector.h"
 #include "openmc/xml_interface.h"
 
@@ -174,8 +174,8 @@ class Mesh {
   //! \param[out] Array of (material index, volume) for desired element
   //! \param[inout] seed Pseudorandom number seed
   //! \return Number of materials within element
-  int material_volumes(int n_sample, int bin, gsl::span<MaterialVolume> volumes,
-    uint64_t* seed) const;
+  int material_volumes(
+    int n_sample, int bin, span<MaterialVolume> volumes, uint64_t* seed) const;
 
   //! Determine volume of materials within a single mesh elemenet
   //
diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h
index 1cc9d297287..329c776d032 100644
--- a/include/openmc/nuclide.h
+++ b/include/openmc/nuclide.h
@@ -7,7 +7,6 @@
 #include <unordered_map>
 #include <utility> // for pair
 
-#include <gsl/gsl-lite.hpp>
 #include <hdf5.h>
 
 #include "openmc/array.h"
@@ -17,6 +16,7 @@
 #include "openmc/particle.h"
 #include "openmc/reaction.h"
 #include "openmc/reaction_product.h"
+#include "openmc/span.h"
 #include "openmc/urr.h"
 #include "openmc/vector.h"
 #include "openmc/wmp.h"
@@ -81,8 +81,8 @@ class Nuclide {
   //! \param[in] energy Energy group boundaries in [eV]
   //! \param[in] flux Flux in each energy group (not normalized per eV)
   //! \return Reaction rate
-  double collapse_rate(int MT, double temperature,
-    gsl::span<const double> energy, gsl::span<const double> flux) const;
+  double collapse_rate(int MT, double temperature, span<const double> energy,
+    span<const double> flux) const;
 
   //============================================================================
   // Data members
@@ -91,7 +91,7 @@ class Nuclide {
   int A_;            //!< Mass number
   int metastable_;   //!< Metastable state
   double awr_;       //!< Atomic weight ratio
-  gsl::index index_; //!< Index in the nuclides array
+  int64_t index_;    //!< Index in the nuclides array
 
   // Temperature dependent cross section data
   vector<double> kTs_;                //!< temperatures in eV (k*T)
@@ -138,7 +138,7 @@ class Nuclide {
   //
   //! \param[in] T Temperature in [K]
   //! \return Temperature index and interpolation factor
-  std::pair<gsl::index, double> find_temperature(double T) const;
+  std::pair<int64_t, double> find_temperature(double T) const;
 
   static int XS_TOTAL;
   static int XS_ABSORPTION;
diff --git a/include/openmc/photon.h b/include/openmc/photon.h
index 9901c8c6db5..1fee6c9f5b6 100644
--- a/include/openmc/photon.h
+++ b/include/openmc/photon.h
@@ -7,7 +7,6 @@
 #include "openmc/vector.h"
 
 #include "xtensor/xtensor.hpp"
-#include <gsl/gsl-lite.hpp>
 #include <hdf5.h>
 
 #include <string>
@@ -61,7 +60,7 @@ class PhotonInteraction {
   // Data members
   std::string name_; //!< Name of element, e.g. "Zr"
   int Z_;            //!< Atomic number
-  gsl::index index_; //!< Index in global elements vector
+  int64_t index_;    //!< Index in global elements vector
 
   // Microscopic cross sections
   xt::xtensor<double, 1> energy_;
diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h
index 93987d9d7ee..d5f04d136d0 100644
--- a/include/openmc/reaction.h
+++ b/include/openmc/reaction.h
@@ -7,10 +7,10 @@
 #include <string>
 
 #include "hdf5.h"
-#include <gsl/gsl-lite.hpp>
 
 #include "openmc/particle_data.h"
 #include "openmc/reaction_product.h"
+#include "openmc/span.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -33,7 +33,7 @@ class Reaction {
   //! \param[in] i_temp Temperature index
   //! \param[in] i_grid Energy grid index
   //! \param[in] interp_factor Interpolation factor between grid points
-  double xs(gsl::index i_temp, gsl::index i_grid, double interp_factor) const;
+  double xs(int64_t i_temp, int64_t i_grid, double interp_factor) const;
 
   //! Calculate cross section
   //
@@ -47,8 +47,8 @@ class Reaction {
   //! \param[in] flux Flux in each energy group (not normalized per eV)
   //! \param[in] grid Nuclide energy grid
   //! \return Reaction rate
-  double collapse_rate(gsl::index i_temp, gsl::span<const double> energy,
-    gsl::span<const double> flux, const vector<double>& grid) const;
+  double collapse_rate(int64_t i_temp, span<const double> energy,
+    span<const double> flux, const vector<double>& grid) const;
 
   //! Cross section at a single temperature
   struct TemperatureXS {
diff --git a/include/openmc/span.h b/include/openmc/span.h
new file mode 100644
index 00000000000..723bccd76e2
--- /dev/null
+++ b/include/openmc/span.h
@@ -0,0 +1,237 @@
+#ifndef OPENMC_SPAN_H
+#define OPENMC_SPAN_H
+#include <cstddef>   // for std::size_t, std::ptrdiff_t
+#include <iterator>  // for std::begin, std::end
+#include <stdexcept> // for std::out_of_range
+#include <type_traits>
+
+#include "openmc/vector.h"
+
+namespace openmc {
+
+template<typename T>
+class span {
+public:
+  using value_type = T;
+  using pointer = T*;
+  using const_pointer = const T*;
+  using reference = T&;
+  using const_reference = const T&;
+  using iterator = T*;
+  using const_iterator = const T*;
+  using size_type = std::size_t;
+  using difference_type = std::ptrdiff_t;
+
+  /**
+   * @brief Default constructor for an empty span.
+   */
+  span() noexcept : data_(nullptr), size_(0) {}
+
+  /**
+   * @brief Constructs a span from a pointer and size.
+   * @param ptr Pointer to the first element.
+   * @param count Number of elements in the span.
+   */
+  span(pointer ptr, size_type count) : data_(ptr), size_(count) {}
+
+  /**
+   * @brief Constructs a span from two pointers marking the span range.
+   * @param first Pointer to the first element.
+   * @param last Pointer past the last element.
+   * @throws std::out_of_range if last < first.
+   */
+  span(pointer first, pointer last) : data_(first), size_(last - first)
+  {
+    if (last < first) {
+      throw std::out_of_range("span: last pointer is before first pointer");
+    }
+  }
+
+  /**
+   * @brief Constructs a span from a non-const std::vector.
+   * @param vec Reference to the vector to create a span from.
+   */
+  span(std::vector<T>& vec) : data_(vec.data()), size_(vec.size()) {}
+
+  /**
+   * @brief Constructs a span from a const std::vector.
+   *
+   * This is handling the semantics that a span<const double> is used
+   * for read-only access into a vector.
+   * @param vec Reference to the const vector to create a span from.
+   */
+  template<typename U = T,
+    typename = std::enable_if_t<std::is_same<U, const T>::value>>
+  span(const std::vector<std::remove_const_t<U>>& vec)
+    : data_(vec.data()), size_(vec.size())
+  {}
+
+  /**
+   * @brief Constructs a read-only span<const T> from a non-const span<T>.
+   */
+  template<typename U = T, typename = std::enable_if_t<std::is_const<U>::value>>
+  span(const span<std::remove_const_t<U>>& other) noexcept
+    : data_(other.data()), size_(other.size())
+  {}
+
+  /**
+   * @brief Access an element without bounds checking.
+   * @param index Index of the element to access.
+   * @return Reference to the accessed element.
+   */
+  reference operator[](size_type index) { return data_[index]; }
+
+  /**
+   * @brief Access an element without bounds checking (const version).
+   * @param index Index of the element to access.
+   * @return Const reference to the accessed element.
+   */
+  const_reference operator[](size_type index) const { return data_[index]; }
+
+  /**
+   * @brief Access an element with bounds checking.
+   * @param index Index of the element to access.
+   * @return Reference to the accessed element.
+   * @throws std::out_of_range if index is out of range.
+   */
+  reference at(size_type index)
+  {
+    if (index >= size_) {
+      throw std::out_of_range("span: index out of range");
+    }
+    return data_[index];
+  }
+
+  /**
+   * @brief Access an element with bounds checking (const version).
+   * @param index Index of the element to access.
+   * @return Const reference to the accessed element.
+   * @throws std::out_of_range if index is out of range.
+   */
+  const_reference at(size_type index) const
+  {
+    if (index >= size_) {
+      throw std::out_of_range("span: index out of range");
+    }
+    return data_[index];
+  }
+
+  /**
+   * @brief Get a pointer to the underlying data.
+   * @return Pointer to the data, or nullptr if the span is empty.
+   */
+  pointer data() noexcept { return data_; }
+
+  /**
+   * @brief Get a const pointer to the underlying data.
+   * @return Const pointer to the data, or nullptr if the span is empty.
+   */
+  const_pointer data() const noexcept { return data_; }
+
+  /**
+   * @brief Get the number of elements in the span.
+   * @return The size of the span.
+   */
+  size_type size() const noexcept { return size_; }
+
+  /**
+   * @brief Check if the span is empty.
+   * @return True if the span is empty, false otherwise.
+   */
+  bool empty() const noexcept { return size_ == 0; }
+
+  /**
+   * @brief Get an iterator to the beginning of the span.
+   * @return Iterator pointing to the first element.
+   */
+  iterator begin() noexcept { return data_; }
+
+  /**
+   * @brief Get a const iterator to the beginning of the span.
+   * @return Const iterator pointing to the first element.
+   */
+  const_iterator begin() const noexcept { return data_; }
+
+  /**
+   * @brief Get a const iterator to the beginning of the span.
+   * @return Const iterator pointing to the first element.
+   */
+  const_iterator cbegin() const noexcept { return data_; }
+
+  /**
+   * @brief Get an iterator to the end of the span.
+   * @return Iterator pointing past the last element.
+   */
+  iterator end() noexcept { return data_ + size_; }
+
+  /**
+   * @brief Get a const iterator to the end of the span.
+   * @return Const iterator pointing past the last element.
+   */
+  const_iterator end() const noexcept { return data_ + size_; }
+
+  /**
+   * @brief Get a const iterator to the end of the span.
+   * @return Const iterator pointing past the last element.
+   */
+  const_iterator cend() const noexcept { return data_ + size_; }
+
+  /**
+   * @brief Access the first element.
+   * @return Reference to the first element.
+   * @throws std::out_of_range if the span is empty.
+   */
+  reference front()
+  {
+    if (empty()) {
+      throw std::out_of_range("span::front(): span is empty");
+    }
+    return data_[0];
+  }
+
+  /**
+   * @brief Access the first element (const version).
+   * @return Const reference to the first element.
+   * @throws std::out_of_range if the span is empty.
+   */
+  const_reference front() const
+  {
+    if (empty()) {
+      throw std::out_of_range("span::front(): span is empty");
+    }
+    return data_[0];
+  }
+
+  /**
+   * @brief Access the last element.
+   * @return Reference to the last element.
+   * @throws std::out_of_range if the span is empty.
+   */
+  reference back()
+  {
+    if (empty()) {
+      throw std::out_of_range("span::back(): span is empty");
+    }
+    return data_[size_ - 1];
+  }
+
+  /**
+   * @brief Access the last element (const version).
+   * @return Const reference to the last element.
+   * @throws std::out_of_range if the span is empty.
+   */
+  const_reference back() const
+  {
+    if (empty()) {
+      throw std::out_of_range("span::back(): span is empty");
+    }
+    return data_[size_ - 1];
+  }
+
+private:
+  pointer data_;
+  size_type size_;
+};
+
+} // namespace openmc
+#endif // OPENMC_SPAN_H
diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h
index f0d41e1697f..d8a3bb0e148 100644
--- a/include/openmc/state_point.h
+++ b/include/openmc/state_point.h
@@ -4,13 +4,12 @@
 #include <cstdint>
 #include <string>
 
-#include <gsl/gsl-lite.hpp>
-
 #include "hdf5.h"
 
 #include "openmc/capi.h"
 #include "openmc/particle.h"
 #include "openmc/shared_array.h"
+#include "openmc/span.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -35,14 +34,14 @@ void load_state_point();
 // can be created by calling calculate_parallel_index_vector on
 // source_bank.size() if such a vector is not already available.
 void write_h5_source_point(const char* filename,
-  gsl::span<SourceSite> source_bank, const vector<int64_t>& bank_index);
+  span<const SourceSite> source_bank, const vector<int64_t>& bank_index);
 
-void write_source_point(std::string, gsl::span<SourceSite> source_bank,
+void write_source_point(std::string, span<const SourceSite> source_bank,
   const vector<int64_t>& bank_index, bool use_mcpl);
 
 // This appends a source bank specification to an HDF5 file
 // that's already open. It is used internally by write_source_point.
-void write_source_bank(hid_t group_id, gsl::span<SourceSite> source_bank,
+void write_source_bank(hid_t group_id, span<const SourceSite> source_bank,
   const vector<int64_t>& bank_index);
 
 void read_source_bank(
diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h
index 210ab284ba9..31dd609ed10 100644
--- a/include/openmc/tallies/filter.h
+++ b/include/openmc/tallies/filter.h
@@ -6,7 +6,6 @@
 #include <unordered_map>
 
 #include "pugixml.hpp"
-#include <gsl/gsl-lite.hpp>
 
 #include "openmc/constants.h"
 #include "openmc/hdf5_interface.h"
@@ -130,7 +129,7 @@ class Filter {
   //! \return Number of bins
   int n_bins() const { return n_bins_; }
 
-  gsl::index index() const { return index_; }
+  int64_t index() const { return index_; }
 
   //----------------------------------------------------------------------------
   // Data members
@@ -140,7 +139,7 @@ class Filter {
 
 private:
   int32_t id_ {C_NONE};
-  gsl::index index_;
+  int64_t index_;
 };
 
 //==============================================================================
diff --git a/include/openmc/tallies/filter_azimuthal.h b/include/openmc/tallies/filter_azimuthal.h
index 37e1ef07356..4853c554590 100644
--- a/include/openmc/tallies/filter_azimuthal.h
+++ b/include/openmc/tallies/filter_azimuthal.h
@@ -4,8 +4,7 @@
 #include "openmc/vector.h"
 #include <string>
 
-#include <gsl/gsl-lite.hpp>
-
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 
 namespace openmc {
@@ -39,7 +38,7 @@ class AzimuthalFilter : public Filter {
   //----------------------------------------------------------------------------
   // Accessors
 
-  void set_bins(gsl::span<double> bins);
+  void set_bins(span<double> bins);
 
 private:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/filter_cell.h b/include/openmc/tallies/filter_cell.h
index b7ea01ce6cb..ac6539466f8 100644
--- a/include/openmc/tallies/filter_cell.h
+++ b/include/openmc/tallies/filter_cell.h
@@ -4,8 +4,7 @@
 #include <cstdint>
 #include <unordered_map>
 
-#include <gsl/gsl-lite.hpp>
-
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -42,7 +41,7 @@ class CellFilter : public Filter {
 
   const vector<int32_t>& cells() const { return cells_; }
 
-  void set_cells(gsl::span<int32_t> cells);
+  void set_cells(span<int32_t> cells);
 
 protected:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/filter_cell_instance.h b/include/openmc/tallies/filter_cell_instance.h
index f500f48896f..69499765b22 100644
--- a/include/openmc/tallies/filter_cell_instance.h
+++ b/include/openmc/tallies/filter_cell_instance.h
@@ -4,9 +4,8 @@
 #include <cstdint>
 #include <unordered_map>
 
-#include <gsl/gsl-lite.hpp>
-
 #include "openmc/cell.h"
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -22,7 +21,7 @@ class CellInstanceFilter : public Filter {
   // Constructors, destructors
 
   CellInstanceFilter() = default;
-  CellInstanceFilter(gsl::span<CellInstance> instances);
+  CellInstanceFilter(span<CellInstance> instances);
   ~CellInstanceFilter() = default;
 
   //----------------------------------------------------------------------------
@@ -47,7 +46,7 @@ class CellInstanceFilter : public Filter {
 
   const std::unordered_set<int32_t>& cells() const { return cells_; }
 
-  void set_cell_instances(gsl::span<CellInstance> instances);
+  void set_cell_instances(span<CellInstance> instances);
 
 private:
   //----------------------------------------------------------------------------
@@ -60,7 +59,7 @@ class CellInstanceFilter : public Filter {
   std::unordered_set<int32_t> cells_;
 
   //! A map from cell/instance indices to filter bin indices.
-  std::unordered_map<CellInstance, gsl::index, CellInstanceHash> map_;
+  std::unordered_map<CellInstance, int64_t, CellInstanceHash> map_;
 
   //! Indicates if filter uses only material-filled cells
   bool material_cells_only_;
diff --git a/include/openmc/tallies/filter_collision.h b/include/openmc/tallies/filter_collision.h
index d2dab70ca06..7d42a5ddd35 100644
--- a/include/openmc/tallies/filter_collision.h
+++ b/include/openmc/tallies/filter_collision.h
@@ -1,9 +1,9 @@
 #ifndef OPENMC_TALLIES_FILTER_COLLISIONS_H
 #define OPENMC_TALLIES_FILTER_COLLISIONS_H
 
-#include <gsl/gsl-lite.hpp>
 #include <unordered_map>
 
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -39,7 +39,7 @@ class CollisionFilter : public Filter {
   // Accessors
 
   const vector<int>& bins() const { return bins_; }
-  void set_bins(gsl::span<const int> bins);
+  void set_bins(span<const int> bins);
 
 protected:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/filter_delayedgroup.h b/include/openmc/tallies/filter_delayedgroup.h
index 71919b2ece5..7d11447ab13 100644
--- a/include/openmc/tallies/filter_delayedgroup.h
+++ b/include/openmc/tallies/filter_delayedgroup.h
@@ -1,8 +1,7 @@
 #ifndef OPENMC_TALLIES_FILTER_DELAYEDGROUP_H
 #define OPENMC_TALLIES_FILTER_DELAYEDGROUP_H
 
-#include <gsl/gsl-lite.hpp>
-
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -42,7 +41,7 @@ class DelayedGroupFilter : public Filter {
 
   const vector<int>& groups() const { return groups_; }
 
-  void set_groups(gsl::span<int> groups);
+  void set_groups(span<int> groups);
 
 private:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/filter_energy.h b/include/openmc/tallies/filter_energy.h
index e35e01a6dae..cf8a8aa0f58 100644
--- a/include/openmc/tallies/filter_energy.h
+++ b/include/openmc/tallies/filter_energy.h
@@ -1,8 +1,7 @@
 #ifndef OPENMC_TALLIES_FILTER_ENERGY_H
 #define OPENMC_TALLIES_FILTER_ENERGY_H
 
-#include <gsl/gsl-lite.hpp>
-
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -38,7 +37,7 @@ class EnergyFilter : public Filter {
   // Accessors
 
   const vector<double>& bins() const { return bins_; }
-  void set_bins(gsl::span<const double> bins);
+  void set_bins(span<const double> bins);
 
   bool matches_transport_groups() const { return matches_transport_groups_; }
 
diff --git a/include/openmc/tallies/filter_energyfunc.h b/include/openmc/tallies/filter_energyfunc.h
index d77ef0fa839..e03c23dda7e 100644
--- a/include/openmc/tallies/filter_energyfunc.h
+++ b/include/openmc/tallies/filter_energyfunc.h
@@ -2,6 +2,7 @@
 #define OPENMC_TALLIES_FILTER_ENERGYFUNC_H
 
 #include "openmc/constants.h"
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -42,7 +43,7 @@ class EnergyFunctionFilter : public Filter {
   const vector<double>& energy() const { return energy_; }
   const vector<double>& y() const { return y_; }
   Interpolation interpolation() const { return interpolation_; }
-  void set_data(gsl::span<const double> energy, gsl::span<const double> y);
+  void set_data(span<const double> energy, span<const double> y);
   void set_interpolation(const std::string& interpolation);
 
 private:
diff --git a/include/openmc/tallies/filter_material.h b/include/openmc/tallies/filter_material.h
index aa5df5b3ab2..ccfe5403d5e 100644
--- a/include/openmc/tallies/filter_material.h
+++ b/include/openmc/tallies/filter_material.h
@@ -4,8 +4,7 @@
 #include <cstdint>
 #include <unordered_map>
 
-#include <gsl/gsl-lite.hpp>
-
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -44,7 +43,7 @@ class MaterialFilter : public Filter {
 
   const vector<int32_t>& materials() const { return materials_; }
 
-  void set_materials(gsl::span<const int32_t> materials);
+  void set_materials(span<const int32_t> materials);
 
 protected:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/filter_mu.h b/include/openmc/tallies/filter_mu.h
index b0ee40eac9f..d6e7f1798e2 100644
--- a/include/openmc/tallies/filter_mu.h
+++ b/include/openmc/tallies/filter_mu.h
@@ -1,8 +1,7 @@
 #ifndef OPENMC_TALLIES_FILTER_MU_H
 #define OPENMC_TALLIES_FILTER_MU_H
 
-#include <gsl/gsl-lite.hpp>
-
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -38,7 +37,7 @@ class MuFilter : public Filter {
   //----------------------------------------------------------------------------
   // Accessors
 
-  void set_bins(gsl::span<double> bins);
+  void set_bins(span<double> bins);
 
 protected:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/filter_musurface.h b/include/openmc/tallies/filter_musurface.h
index 2ca19e3a259..fa6816836d9 100644
--- a/include/openmc/tallies/filter_musurface.h
+++ b/include/openmc/tallies/filter_musurface.h
@@ -1,8 +1,6 @@
 #ifndef OPENMC_TALLIES_FILTER_MU_SURFACE_H
 #define OPENMC_TALLIES_FILTER_MU_SURFACE_H
 
-#include <gsl/gsl-lite.hpp>
-
 #include "openmc/tallies/filter_mu.h"
 #include "openmc/vector.h"
 
diff --git a/include/openmc/tallies/filter_particle.h b/include/openmc/tallies/filter_particle.h
index a181d5cee64..863a6d282fe 100644
--- a/include/openmc/tallies/filter_particle.h
+++ b/include/openmc/tallies/filter_particle.h
@@ -2,6 +2,7 @@
 #define OPENMC_TALLIES_FILTER_PARTICLE_H
 
 #include "openmc/particle.h"
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -38,7 +39,7 @@ class ParticleFilter : public Filter {
 
   const vector<ParticleType>& particles() const { return particles_; }
 
-  void set_particles(gsl::span<ParticleType> particles);
+  void set_particles(span<ParticleType> particles);
 
 private:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/filter_polar.h b/include/openmc/tallies/filter_polar.h
index 78bb25aa45e..c7c73c89f54 100644
--- a/include/openmc/tallies/filter_polar.h
+++ b/include/openmc/tallies/filter_polar.h
@@ -3,8 +3,7 @@
 
 #include <cmath>
 
-#include <gsl/gsl-lite.hpp>
-
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -39,7 +38,7 @@ class PolarFilter : public Filter {
   //----------------------------------------------------------------------------
   // Accessors
 
-  void set_bins(gsl::span<double> bins);
+  void set_bins(span<double> bins);
 
 private:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/filter_sph_harm.h b/include/openmc/tallies/filter_sph_harm.h
index 5d4a4bd999c..a6904c30103 100644
--- a/include/openmc/tallies/filter_sph_harm.h
+++ b/include/openmc/tallies/filter_sph_harm.h
@@ -3,8 +3,6 @@
 
 #include <string>
 
-#include <gsl/gsl-lite.hpp>
-
 #include "openmc/tallies/filter.h"
 
 namespace openmc {
@@ -46,7 +44,7 @@ class SphericalHarmonicsFilter : public Filter {
 
   SphericalHarmonicsCosine cosine() const { return cosine_; }
 
-  void set_cosine(gsl::cstring_span cosine);
+  void set_cosine(const std::string& cosine);
 
 private:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/filter_surface.h b/include/openmc/tallies/filter_surface.h
index 3537f1cc748..e78243f5ff2 100644
--- a/include/openmc/tallies/filter_surface.h
+++ b/include/openmc/tallies/filter_surface.h
@@ -4,8 +4,7 @@
 #include <cstdint>
 #include <unordered_map>
 
-#include <gsl/gsl-lite.hpp>
-
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -40,7 +39,7 @@ class SurfaceFilter : public Filter {
   //----------------------------------------------------------------------------
   // Accessors
 
-  void set_surfaces(gsl::span<int32_t> surfaces);
+  void set_surfaces(span<int32_t> surfaces);
 
 private:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/filter_time.h b/include/openmc/tallies/filter_time.h
index 105ef9880a7..3ce557abda6 100644
--- a/include/openmc/tallies/filter_time.h
+++ b/include/openmc/tallies/filter_time.h
@@ -1,8 +1,7 @@
 #ifndef OPENMC_TALLIES_FILTER_TIME_H
 #define OPENMC_TALLIES_FILTER_TIME_H
 
-#include <gsl/gsl-lite.hpp>
-
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -38,7 +37,7 @@ class TimeFilter : public Filter {
   // Accessors
 
   const vector<double>& bins() const { return bins_; }
-  void set_bins(gsl::span<const double> bins);
+  void set_bins(span<const double> bins);
 
 protected:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/filter_universe.h b/include/openmc/tallies/filter_universe.h
index d4894353bc9..461434ec441 100644
--- a/include/openmc/tallies/filter_universe.h
+++ b/include/openmc/tallies/filter_universe.h
@@ -4,8 +4,7 @@
 #include <cstdint>
 #include <unordered_map>
 
-#include <gsl/gsl-lite.hpp>
-
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -40,7 +39,7 @@ class UniverseFilter : public Filter {
   //----------------------------------------------------------------------------
   // Accessors
 
-  void set_universes(gsl::span<int32_t> universes);
+  void set_universes(span<int32_t> universes);
 
 private:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
index 2f0cddcf0ae..48f678ced0f 100644
--- a/include/openmc/tallies/tally.h
+++ b/include/openmc/tallies/tally.h
@@ -3,6 +3,7 @@
 
 #include "openmc/constants.h"
 #include "openmc/memory.h" // for unique_ptr
+#include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/trigger.h"
 #include "openmc/vector.h"
@@ -10,7 +11,6 @@
 #include "pugixml.hpp"
 #include "xtensor/xfixed.hpp"
 #include "xtensor/xtensor.hpp"
-#include <gsl/gsl-lite.hpp>
 
 #include <string>
 #include <unordered_map>
@@ -93,7 +93,7 @@ class Tally {
   //! \brief Check if this tally has a specified type of filter
   bool has_filter(FilterType filter_type) const;
 
-  void set_filters(gsl::span<Filter*> filters);
+  void set_filters(span<Filter*> filters);
 
   //! Given already-set filters, set the stride lengths
   void set_strides();
@@ -192,7 +192,7 @@ class Tally {
   //! Whether to multiply by atom density for reaction rates
   bool multiply_density_ {true};
 
-  gsl::index index_;
+  int64_t index_;
 };
 
 //==============================================================================
diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h
index 376b8c707dd..fa8d3d65ece 100644
--- a/include/openmc/volume_calc.h
+++ b/include/openmc/volume_calc.h
@@ -14,7 +14,6 @@
 
 #include "pugixml.hpp"
 #include "xtensor/xtensor.hpp"
-#include <gsl/gsl-lite.hpp>
 #ifdef _OPENMP
 #include <omp.h>
 #endif
diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h
index 9852fe679d0..1ec8cc512df 100644
--- a/include/openmc/weight_windows.h
+++ b/include/openmc/weight_windows.h
@@ -4,7 +4,6 @@
 #include <cstdint>
 #include <unordered_map>
 
-#include <gsl/gsl-lite.hpp>
 #include <hdf5.h>
 #include <pugixml.hpp>
 
@@ -12,6 +11,7 @@
 #include "openmc/memory.h"
 #include "openmc/mesh.h"
 #include "openmc/particle.h"
+#include "openmc/span.h"
 #include "openmc/tallies/tally.h"
 #include "openmc/vector.h"
 
@@ -106,7 +106,7 @@ class WeightWindows {
   //! Set the weight window ID
   void set_id(int32_t id = -1);
 
-  void set_energy_bounds(gsl::span<const double> bounds);
+  void set_energy_bounds(span<const double> bounds);
 
   void set_mesh(const std::unique_ptr<Mesh>& mesh);
 
@@ -149,9 +149,9 @@ class WeightWindows {
   void set_bounds(const xt::xtensor<double, 2>& lower_bounds, double ratio);
 
   void set_bounds(
-    gsl::span<const double> lower_bounds, gsl::span<const double> upper_bounds);
+    span<const double> lower_bounds, span<const double> upper_bounds);
 
-  void set_bounds(gsl::span<const double> lower_bounds, double ratio);
+  void set_bounds(span<const double> lower_bounds, double ratio);
 
   void set_particle_type(ParticleType p_type);
 
@@ -193,8 +193,8 @@ class WeightWindows {
 private:
   //----------------------------------------------------------------------------
   // Data members
-  int32_t id_;       //!< Unique ID
-  gsl::index index_; //!< Index into weight windows vector
+  int32_t id_;    //!< Unique ID
+  int64_t index_; //!< Index into weight windows vector
   ParticleType particle_type_ {
     ParticleType::neutron};      //!< Particle type to apply weight windows to
   vector<double> energy_bounds_; //!< Energy boundaries [eV]
diff --git a/src/cell.cpp b/src/cell.cpp
index 88876678706..a2170aed11f 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -2,6 +2,7 @@
 #include "openmc/cell.h"
 
 #include <algorithm>
+#include <cassert>
 #include <cctype>
 #include <cmath>
 #include <iterator>
@@ -10,7 +11,6 @@
 #include <string>
 
 #include <fmt/core.h>
-#include <gsl/gsl-lite.hpp>
 
 #include "openmc/capi.h"
 #include "openmc/constants.h"
@@ -137,7 +137,7 @@ void Cell::set_temperature(double T, int32_t instance, bool set_contained)
     auto contained_cells = this->get_contained_cells(instance);
     for (const auto& entry : contained_cells) {
       auto& cell = model::cells[entry.first];
-      Expects(cell->type_ == Fill::MATERIAL);
+      assert(cell->type_ == Fill::MATERIAL);
       auto& instances = entry.second;
       for (auto instance : instances) {
         cell->set_temperature(T, instance);
@@ -179,7 +179,7 @@ void Cell::import_properties_hdf5(hid_t group)
   // Modify temperatures for the cell
   sqrtkT_.clear();
   sqrtkT_.resize(temps.size());
-  for (gsl::index i = 0; i < temps.size(); ++i) {
+  for (int64_t i = 0; i < temps.size(); ++i) {
     this->set_temperature(temps[i], i);
   }
 
@@ -578,7 +578,7 @@ void Region::apply_demorgan(
 //! precedence than unions using parentheses.
 //==============================================================================
 
-gsl::index Region::add_parentheses(gsl::index start)
+int64_t Region::add_parentheses(int64_t start)
 {
   int32_t start_token = expression_[start];
   // Add left parenthesis and set new position to be after parenthesis
@@ -650,7 +650,7 @@ void Region::add_precedence()
   int32_t current_op = 0;
   std::size_t current_dist = 0;
 
-  for (gsl::index i = 0; i < expression_.size(); i++) {
+  for (int64_t i = 0; i < expression_.size(); i++) {
     int32_t token = expression_[i];
 
     if (token == OP_UNION || token == OP_INTERSECTION) {
@@ -946,7 +946,7 @@ BoundingBox Region::bounding_box_complex(vector<int32_t> postfix) const
     }
   }
 
-  Ensures(i_stack == 0);
+  assert(i_stack == 0);
   return stack.front();
 }
 
@@ -1218,8 +1218,8 @@ struct ParentCell {
              lattice_index < other.lattice_index);
   }
 
-  gsl::index cell_index;
-  gsl::index lattice_index;
+  int64_t cell_index;
+  int64_t lattice_index;
 };
 
 //! Structure used to insert ParentCell into hashed STL data structures
diff --git a/src/dagmc.cpp b/src/dagmc.cpp
index b79676c3626..5be156bc903 100644
--- a/src/dagmc.cpp
+++ b/src/dagmc.cpp
@@ -1,5 +1,7 @@
 #include "openmc/dagmc.h"
 
+#include <cassert>
+
 #include "openmc/constants.h"
 #include "openmc/container_util.h"
 #include "openmc/error.h"
@@ -760,7 +762,7 @@ Direction DAGSurface::normal(Position r) const
 
 Direction DAGSurface::reflect(Position r, Direction u, GeometryState* p) const
 {
-  Expects(p);
+  assert(p);
   p->history().reset_to_last_intersection();
   moab::ErrorCode rval;
   moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
diff --git a/src/distribution.cpp b/src/distribution.cpp
index a6b4acd58b1..cfc9d311038 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -8,8 +8,6 @@
 #include <stdexcept> // for runtime_error
 #include <string>    // for string, stod
 
-#include <gsl/gsl-lite.hpp>
-
 #include "openmc/error.h"
 #include "openmc/math_functions.h"
 #include "openmc/random_dist.h"
@@ -30,12 +28,12 @@ DiscreteIndex::DiscreteIndex(pugi::xml_node node)
   assign({params.data() + n, n});
 }
 
-DiscreteIndex::DiscreteIndex(gsl::span<const double> p)
+DiscreteIndex::DiscreteIndex(span<const double> p)
 {
   assign(p);
 }
 
-void DiscreteIndex::assign(gsl::span<const double> p)
+void DiscreteIndex::assign(span<const double> p)
 {
   prob_.assign(p.begin(), p.end());
 
@@ -413,7 +411,7 @@ double Mixture::sample(uint64_t* seed) const
     p, [](const DistPair& pair, double p) { return pair.first < p; });
 
   // This should not happen. Catch it
-  Ensures(it != distribution_.cend());
+  assert(it != distribution_.cend());
 
   // Sample the chosen distribution
   return it->second->sample(seed);
diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp
index 5c193a95d29..d2b0f413bd1 100644
--- a/src/distribution_spatial.cpp
+++ b/src/distribution_spatial.cpp
@@ -262,7 +262,7 @@ MeshSpatial::MeshSpatial(pugi::xml_node node)
   elem_idx_dist_.assign(strengths);
 }
 
-MeshSpatial::MeshSpatial(int32_t mesh_idx, gsl::span<const double> strengths)
+MeshSpatial::MeshSpatial(int32_t mesh_idx, span<const double> strengths)
   : mesh_idx_(mesh_idx)
 {
   check_element_types();
@@ -331,7 +331,7 @@ PointCloud::PointCloud(pugi::xml_node node)
 }
 
 PointCloud::PointCloud(
-  std::vector<Position> point_cloud, gsl::span<const double> strengths)
+  std::vector<Position> point_cloud, span<const double> strengths)
 {
   point_cloud_.assign(point_cloud.begin(), point_cloud.end());
   point_idx_dist_.assign(strengths);
diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp
index 050d4db968c..2f8a5574337 100644
--- a/src/geometry_aux.cpp
+++ b/src/geometry_aux.cpp
@@ -344,7 +344,7 @@ void prepare_distribcell(const std::vector<int32_t>* user_distribcells)
 
   // By default, add material cells to the list of distributed cells
   if (settings::material_cell_offsets) {
-    for (gsl::index i = 0; i < model::cells.size(); ++i) {
+    for (int64_t i = 0; i < model::cells.size(); ++i) {
       if (model::cells[i]->type_ == Fill::MATERIAL)
         distribcells.insert(i);
     }
diff --git a/src/material.cpp b/src/material.cpp
index aad7008e703..3ba9ab96c00 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -1,6 +1,7 @@
 #include "openmc/material.h"
 
 #include <algorithm> // for min, max, sort, fill
+#include <cassert>
 #include <cmath>
 #include <iterator>
 #include <sstream>
@@ -933,7 +934,7 @@ void Material::calculate_photon_xs(Particle& p) const
 
 void Material::set_id(int32_t id)
 {
-  Expects(id >= 0 || id == C_NONE);
+  assert(id >= 0 || id == C_NONE);
 
   // Clear entry in material map if an ID was already assigned before
   if (id_ != C_NONE) {
@@ -961,9 +962,9 @@ void Material::set_id(int32_t id)
   model::material_map[id] = index_;
 }
 
-void Material::set_density(double density, gsl::cstring_span units)
+void Material::set_density(double density, const std::string& units)
 {
-  Expects(density >= 0.0);
+  assert(density >= 0.0);
 
   if (nuclide_.empty()) {
     throw std::runtime_error {"No nuclides exist in material yet."};
@@ -1006,8 +1007,8 @@ void Material::set_densities(
   const vector<std::string>& name, const vector<double>& density)
 {
   auto n = name.size();
-  Expects(n > 0);
-  Expects(n == density.size());
+  assert(n > 0);
+  assert(n == density.size());
 
   if (n != nuclide_.size()) {
     nuclide_.resize(n);
@@ -1017,7 +1018,7 @@ void Material::set_densities(
   }
 
   double sum_density = 0.0;
-  for (gsl::index i = 0; i < n; ++i) {
+  for (int64_t i = 0; i < n; ++i) {
     const auto& nuc {name[i]};
     if (data::nuclide_map.find(nuc) == data::nuclide_map.end()) {
       int err = openmc_load_nuclide(nuc.c_str(), nullptr, 0);
@@ -1026,7 +1027,7 @@ void Material::set_densities(
     }
 
     nuclide_[i] = data::nuclide_map.at(nuc);
-    Expects(density[i] > 0.0);
+    assert(density[i] > 0.0);
     atom_density_(i) = density[i];
     sum_density += density[i];
 
diff --git a/src/mcpl_interface.cpp b/src/mcpl_interface.cpp
index 5c3df026ce5..8f52861c9a7 100644
--- a/src/mcpl_interface.cpp
+++ b/src/mcpl_interface.cpp
@@ -112,7 +112,7 @@ vector<SourceSite> mcpl_source_sites(std::string path)
 
 #ifdef OPENMC_MCPL
 void write_mcpl_source_bank(mcpl_outfile_t file_id,
-  gsl::span<SourceSite> source_bank, const vector<int64_t>& bank_index)
+  span<const SourceSite> source_bank, const vector<int64_t>& bank_index)
 {
   int64_t dims_size = settings::n_particles;
   int64_t count_size = simulation::work_per_rank;
@@ -189,7 +189,7 @@ void write_mcpl_source_bank(mcpl_outfile_t file_id,
 //==============================================================================
 
 void write_mcpl_source_point(const char* filename,
-  gsl::span<SourceSite> source_bank, const vector<int64_t>& bank_index)
+  span<const SourceSite> source_bank, const vector<int64_t>& bank_index)
 {
   std::string filename_(filename);
   const auto extension = get_file_extension(filename_);
diff --git a/src/mesh.cpp b/src/mesh.cpp
index c572d255a36..9ac60c8644e 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -1,8 +1,8 @@
 #include "openmc/mesh.h"
 #include <algorithm> // for copy, equal, min, min_element
-#include <cmath>     // for ceil
-#include <cstddef>   // for size_t
-#include <gsl/gsl-lite.hpp>
+#include <cassert>
+#include <cmath>   // for ceil
+#include <cstddef> // for size_t
 #include <string>
 
 #ifdef OPENMC_MPI
@@ -113,7 +113,7 @@ Mesh::Mesh(pugi::xml_node node)
 
 void Mesh::set_id(int32_t id)
 {
-  Expects(id >= 0 || id == C_NONE);
+  assert(id >= 0 || id == C_NONE);
 
   // Clear entry in mesh map in case one was already assigned
   if (id_ != C_NONE) {
@@ -151,7 +151,7 @@ vector<double> Mesh::volumes() const
 }
 
 int Mesh::material_volumes(
-  int n_sample, int bin, gsl::span<MaterialVolume> result, uint64_t* seed) const
+  int n_sample, int bin, span<MaterialVolume> result, uint64_t* seed) const
 {
   vector<int32_t> materials;
   vector<int64_t> hits;
@@ -3016,7 +3016,7 @@ void LibMesh::initialize()
     bin_to_elem_map_.reserve(m_->n_active_elem());
     elem_to_bin_map_.resize(m_->n_elem(), -1);
     for (auto it = m_->active_elements_begin(); it != m_->active_elements_end();
-         it++) {
+      it++) {
       auto elem = *it;
 
       bin_to_elem_map_.push_back(elem->id());
@@ -3175,7 +3175,7 @@ void LibMesh::set_score_data(const std::string& var_name,
   unsigned int std_dev_num = variable_map_.at(std_dev_name);
 
   for (auto it = m_->local_elements_begin(); it != m_->local_elements_end();
-       it++) {
+    it++) {
     if (!(*it)->active()) {
       continue;
     }
@@ -3185,13 +3185,13 @@ void LibMesh::set_score_data(const std::string& var_name,
     // set value
     vector<libMesh::dof_id_type> value_dof_indices;
     dof_map.dof_indices(*it, value_dof_indices, value_num);
-    Ensures(value_dof_indices.size() == 1);
+    assert(value_dof_indices.size() == 1);
     eqn_sys.solution->set(value_dof_indices[0], values.at(bin));
 
     // set std dev
     vector<libMesh::dof_id_type> std_dev_dof_indices;
     dof_map.dof_indices(*it, std_dev_dof_indices, std_dev_num);
-    Ensures(std_dev_dof_indices.size() == 1);
+    assert(std_dev_dof_indices.size() == 1);
     eqn_sys.solution->set(std_dev_dof_indices[0], std_dev.at(bin));
   }
 }
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index f720f848bc2..d78b7a0101d 100644
--- a/src/nuclide.cpp
+++ b/src/nuclide.cpp
@@ -21,7 +21,8 @@
 #include "xtensor/xview.hpp"
 
 #include <algorithm> // for sort, min_element
-#include <string>    // for to_string, stoi
+#include <cassert>
+#include <string> // for to_string, stoi
 
 namespace openmc {
 
@@ -999,19 +1000,19 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
   }
 }
 
-std::pair<gsl::index, double> Nuclide::find_temperature(double T) const
+std::pair<int64_t, double> Nuclide::find_temperature(double T) const
 {
-  Expects(T >= 0.0);
+  assert(T >= 0.0);
 
   // Determine temperature index
-  gsl::index i_temp = 0;
+  int64_t i_temp = 0;
   double f = 0.0;
   double kT = K_BOLTZMANN * T;
-  gsl::index n = kTs_.size();
+  int64_t n = kTs_.size();
   switch (settings::temperature_method) {
   case TemperatureMethod::NEAREST: {
     double max_diff = INFTY;
-    for (gsl::index t = 0; t < n; ++t) {
+    for (int64_t t = 0; t < n; ++t) {
       double diff = std::abs(kTs_[t] - kT);
       if (diff < max_diff) {
         i_temp = t;
@@ -1038,17 +1039,17 @@ std::pair<gsl::index, double> Nuclide::find_temperature(double T) const
     f = (kT - kTs_[i_temp]) / (kTs_[i_temp + 1] - kTs_[i_temp]);
   }
 
-  Ensures(i_temp >= 0 && i_temp < n);
+  assert(i_temp >= 0 && i_temp < n);
 
   return {i_temp, f};
 }
 
 double Nuclide::collapse_rate(int MT, double temperature,
-  gsl::span<const double> energy, gsl::span<const double> flux) const
+  span<const double> energy, span<const double> flux) const
 {
-  Expects(MT > 0);
-  Expects(energy.size() > 0);
-  Expects(energy.size() == flux.size() + 1);
+  assert(MT > 0);
+  assert(energy.size() > 0);
+  assert(energy.size() == flux.size() + 1);
 
   int i_rx = reaction_index_[MT];
   if (i_rx < 0)
@@ -1056,7 +1057,7 @@ double Nuclide::collapse_rate(int MT, double temperature,
   const auto& rx = reactions_[i_rx];
 
   // Determine temperature index
-  gsl::index i_temp;
+  int64_t i_temp;
   double f;
   std::tie(i_temp, f) = this->find_temperature(temperature);
 
diff --git a/src/reaction.cpp b/src/reaction.cpp
index 0643cb08bdd..7edf4ff21a7 100644
--- a/src/reaction.cpp
+++ b/src/reaction.cpp
@@ -65,8 +65,7 @@ Reaction::Reaction(hid_t group, const vector<int>& temperatures)
   }
 }
 
-double Reaction::xs(
-  gsl::index i_temp, gsl::index i_grid, double interp_factor) const
+double Reaction::xs(int64_t i_temp, int64_t i_grid, double interp_factor) const
 {
   // If energy is below threshold, return 0. Otherwise interpolate between
   // nearest grid points
@@ -82,9 +81,8 @@ double Reaction::xs(const NuclideMicroXS& micro) const
   return this->xs(micro.index_temp, micro.index_grid, micro.interp_factor);
 }
 
-double Reaction::collapse_rate(gsl::index i_temp,
-  gsl::span<const double> energy, gsl::span<const double> flux,
-  const vector<double>& grid) const
+double Reaction::collapse_rate(int64_t i_temp, span<const double> energy,
+  span<const double> flux, const vector<double>& grid) const
 {
   // Find index corresponding to first energy
   const auto& xs = xs_[i_temp].value;
diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp
index 0b8e1ab42db..8b9e8737c66 100644
--- a/src/secondary_thermal.cpp
+++ b/src/secondary_thermal.cpp
@@ -4,10 +4,9 @@
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 
-#include <gsl/gsl-lite.hpp>
-
 #include "xtensor/xview.hpp"
 
+#include <cassert>
 #include <cmath> // for log, exp
 
 namespace openmc {
@@ -40,7 +39,7 @@ void CoherentElasticAE::sample(
 
   const auto& energies {xs_.bragg_edges()};
 
-  Expects(E_in >= energies.front());
+  assert(E_in >= energies.front());
 
   const int i = lower_bound_index(energies.begin(), energies.end(), E_in);
 
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 09b3764732b..cf74e05d14a 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -439,7 +439,7 @@ void finalize_batch()
       int w = std::to_string(settings::n_max_batches).size();
       std::string source_point_filename = fmt::format("{0}source.{1:0{2}}",
         settings::path_output, simulation::current_batch, w);
-      gsl::span<SourceSite> bankspan(simulation::source_bank);
+      span<const SourceSite> bankspan(simulation::source_bank);
       write_source_point(source_point_filename, bankspan,
         simulation::work_index, settings::source_mcpl_write);
     }
@@ -448,7 +448,7 @@ void finalize_batch()
     if (settings::source_latest) {
       // note: correct file extension appended automatically
       auto filename = settings::path_output + "source";
-      gsl::span<SourceSite> bankspan(simulation::source_bank);
+      span<const SourceSite> bankspan(simulation::source_bank);
       write_source_point(filename.c_str(), bankspan, simulation::work_index,
         settings::source_mcpl_write);
     }
@@ -470,7 +470,7 @@ void finalize_batch()
       // Get span of source bank and calculate parallel index vector
       auto surf_work_index = mpi::calculate_parallel_index_vector(
         simulation::surf_source_bank.size());
-      gsl::span<SourceSite> surfbankspan(simulation::surf_source_bank.begin(),
+      span<const SourceSite> surfbankspan(simulation::surf_source_bank.begin(),
         simulation::surf_source_bank.size());
 
       // Write surface source file
diff --git a/src/state_point.cpp b/src/state_point.cpp
index adc026fa3c2..0285a902ac1 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -569,8 +569,9 @@ hid_t h5banktype()
   return banktype;
 }
 
-void write_source_point(std::string filename, gsl::span<SourceSite> source_bank,
-  const vector<int64_t>& bank_index, bool use_mcpl)
+void write_source_point(std::string filename,
+  span<const SourceSite> source_bank, const vector<int64_t>& bank_index,
+  bool use_mcpl)
 {
   std::string ext = use_mcpl ? "mcpl" : "h5";
   write_message("Creating source file {}.{} with {} particles ...", filename,
@@ -585,7 +586,7 @@ void write_source_point(std::string filename, gsl::span<SourceSite> source_bank,
 }
 
 void write_h5_source_point(const char* filename,
-  gsl::span<SourceSite> source_bank, const vector<int64_t>& bank_index)
+  span<const SourceSite> source_bank, const vector<int64_t>& bank_index)
 {
   // When using parallel HDF5, the file is written to collectively by all
   // processes. With MPI-only, the file is opened and written by the master
@@ -622,7 +623,7 @@ void write_h5_source_point(const char* filename,
     file_close(file_id);
 }
 
-void write_source_bank(hid_t group_id, gsl::span<SourceSite> source_bank,
+void write_source_bank(hid_t group_id, span<const SourceSite> source_bank,
   const vector<int64_t>& bank_index)
 {
   hid_t banktype = h5banktype();
diff --git a/src/surface.cpp b/src/surface.cpp
index 50ef2a12830..a6bf0cd2b23 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -7,7 +7,6 @@
 #include <utility>
 
 #include <fmt/core.h>
-#include <gsl/gsl-lite.hpp>
 
 #include "openmc/array.h"
 #include "openmc/container_util.h"
diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp
index 074212db44d..0ce1ee20b37 100644
--- a/src/tallies/filter.cpp
+++ b/src/tallies/filter.cpp
@@ -1,7 +1,8 @@
 #include "openmc/tallies/filter.h"
 
 #include <algorithm> // for max
-#include <cstring>   // for strcpy
+#include <cassert>
+#include <cstring> // for strcpy
 #include <string>
 
 #include <fmt/core.h>
@@ -162,7 +163,7 @@ Filter* Filter::create(const std::string& type, int32_t id)
 
 void Filter::set_id(int32_t id)
 {
-  Expects(id >= 0 || id == C_NONE);
+  assert(id >= 0 || id == C_NONE);
 
   // Clear entry in filter map if an ID was already assigned before
   if (id_ != C_NONE) {
diff --git a/src/tallies/filter_azimuthal.cpp b/src/tallies/filter_azimuthal.cpp
index e77aa8bdc1e..6525f326d34 100644
--- a/src/tallies/filter_azimuthal.cpp
+++ b/src/tallies/filter_azimuthal.cpp
@@ -34,14 +34,14 @@ void AzimuthalFilter::from_xml(pugi::xml_node node)
   this->set_bins(bins);
 }
 
-void AzimuthalFilter::set_bins(gsl::span<double> bins)
+void AzimuthalFilter::set_bins(span<double> bins)
 {
   // Clear existing bins
   bins_.clear();
   bins_.reserve(bins.size());
 
   // Copy bins, ensuring they are valid
-  for (gsl::index i = 0; i < bins.size(); ++i) {
+  for (int64_t i = 0; i < bins.size(); ++i) {
     if (i > 0 && bins[i] <= bins[i - 1]) {
       throw std::runtime_error {
         "Azimuthal bins must be monotonically increasing."};
diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp
index 794d2ae08b0..b545801ea8e 100644
--- a/src/tallies/filter_cell.cpp
+++ b/src/tallies/filter_cell.cpp
@@ -1,5 +1,7 @@
 #include "openmc/tallies/filter_cell.h"
 
+#include <cassert>
+
 #include <fmt/core.h>
 
 #include "openmc/capi.h"
@@ -25,7 +27,7 @@ void CellFilter::from_xml(pugi::xml_node node)
   this->set_cells(cells);
 }
 
-void CellFilter::set_cells(gsl::span<int32_t> cells)
+void CellFilter::set_cells(span<int32_t> cells)
 {
   // Clear existing cells
   cells_.clear();
@@ -34,8 +36,8 @@ void CellFilter::set_cells(gsl::span<int32_t> cells)
 
   // Update cells and mapping
   for (auto& index : cells) {
-    Expects(index >= 0);
-    Expects(index < model::cells.size());
+    assert(index >= 0);
+    assert(index < model::cells.size());
     cells_.push_back(index);
     map_[index] = cells_.size() - 1;
   }
diff --git a/src/tallies/filter_cell_instance.cpp b/src/tallies/filter_cell_instance.cpp
index 3e78a5bbed9..0634175dd84 100644
--- a/src/tallies/filter_cell_instance.cpp
+++ b/src/tallies/filter_cell_instance.cpp
@@ -1,5 +1,6 @@
 #include "openmc/tallies/filter_cell_instance.h"
 
+#include <cassert>
 #include <string>
 
 #include <fmt/core.h>
@@ -12,7 +13,7 @@
 
 namespace openmc {
 
-CellInstanceFilter::CellInstanceFilter(gsl::span<CellInstance> instances)
+CellInstanceFilter::CellInstanceFilter(span<CellInstance> instances)
 {
   this->set_cell_instances(instances);
 }
@@ -21,26 +22,26 @@ void CellInstanceFilter::from_xml(pugi::xml_node node)
 {
   // Get cell IDs/instances
   auto cells = get_node_array<int32_t>(node, "bins");
-  Expects(cells.size() % 2 == 0);
+  assert(cells.size() % 2 == 0);
 
   // Convert into vector of CellInstance
   vector<CellInstance> instances;
-  for (gsl::index i = 0; i < cells.size() / 2; ++i) {
+  for (int64_t i = 0; i < cells.size() / 2; ++i) {
     int32_t cell_id = cells[2 * i];
-    gsl::index instance = cells[2 * i + 1];
+    int64_t instance = cells[2 * i + 1];
     auto search = model::cell_map.find(cell_id);
     if (search == model::cell_map.end()) {
       throw std::runtime_error {fmt::format(
         "Could not find cell {} specified on tally filter.", cell_id)};
     }
-    gsl::index index = search->second;
+    int64_t index = search->second;
     instances.push_back({index, instance});
   }
 
   this->set_cell_instances(instances);
 }
 
-void CellInstanceFilter::set_cell_instances(gsl::span<CellInstance> instances)
+void CellInstanceFilter::set_cell_instances(span<CellInstance> instances)
 {
   // Clear existing cells
   cell_instances_.clear();
@@ -50,8 +51,8 @@ void CellInstanceFilter::set_cell_instances(gsl::span<CellInstance> instances)
 
   // Update cells and mapping
   for (auto& x : instances) {
-    Expects(x.index_cell >= 0);
-    Expects(x.index_cell < model::cells.size());
+    assert(x.index_cell >= 0);
+    assert(x.index_cell < model::cells.size());
     cell_instances_.push_back(x);
     cells_.insert(x.index_cell);
     map_[x] = cell_instances_.size() - 1;
@@ -72,8 +73,8 @@ void CellInstanceFilter::set_cell_instances(gsl::span<CellInstance> instances)
 void CellInstanceFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
-  gsl::index index_cell = p.lowest_coord().cell;
-  gsl::index instance = p.cell_instance();
+  int64_t index_cell = p.lowest_coord().cell;
+  int64_t instance = p.cell_instance();
 
   if (cells_.count(index_cell) > 0) {
     auto search = map_.find({index_cell, instance});
@@ -88,13 +89,13 @@ void CellInstanceFilter::get_all_bins(
     return;
 
   for (int i = 0; i < p.n_coord() - 1; i++) {
-    gsl::index index_cell = p.coord(i).cell;
+    int64_t index_cell = p.coord(i).cell;
     // if this cell isn't used on the filter, move on
     if (cells_.count(index_cell) == 0)
       continue;
 
     // if this cell is used in the filter, check the instance as well
-    gsl::index instance = cell_instance_at_level(p, i);
+    int64_t instance = cell_instance_at_level(p, i);
     auto search = map_.find({index_cell, instance});
     if (search != map_.end()) {
       match.bins_.push_back(search->second);
@@ -108,7 +109,7 @@ void CellInstanceFilter::to_statepoint(hid_t filter_group) const
   Filter::to_statepoint(filter_group);
   size_t n = cell_instances_.size();
   xt::xtensor<size_t, 2> data({n, 2});
-  for (gsl::index i = 0; i < n; ++i) {
+  for (int64_t i = 0; i < n; ++i) {
     const auto& x = cell_instances_[i];
     data(i, 0) = model::cells[x.index_cell]->id_;
     data(i, 1) = x.instance;
diff --git a/src/tallies/filter_collision.cpp b/src/tallies/filter_collision.cpp
index fbb186a2381..c614c3c83ce 100644
--- a/src/tallies/filter_collision.cpp
+++ b/src/tallies/filter_collision.cpp
@@ -19,7 +19,7 @@ void CollisionFilter::from_xml(pugi::xml_node node)
   this->set_bins(bins);
 }
 
-void CollisionFilter::set_bins(gsl::span<const int> bins)
+void CollisionFilter::set_bins(span<const int> bins)
 {
   // Clear existing bins
   bins_.clear();
@@ -27,7 +27,7 @@ void CollisionFilter::set_bins(gsl::span<const int> bins)
   map_.clear();
 
   // Copy bins
-  for (gsl::index i = 0; i < bins.size(); ++i) {
+  for (int64_t i = 0; i < bins.size(); ++i) {
     bins_.push_back(bins[i]);
     map_[bins[i]] = i;
   }
diff --git a/src/tallies/filter_delayedgroup.cpp b/src/tallies/filter_delayedgroup.cpp
index c6ec217666a..01e39e554a9 100644
--- a/src/tallies/filter_delayedgroup.cpp
+++ b/src/tallies/filter_delayedgroup.cpp
@@ -11,7 +11,7 @@ void DelayedGroupFilter::from_xml(pugi::xml_node node)
   this->set_groups(groups);
 }
 
-void DelayedGroupFilter::set_groups(gsl::span<int> groups)
+void DelayedGroupFilter::set_groups(span<int> groups)
 {
   // Clear existing groups
   groups_.clear();
diff --git a/src/tallies/filter_distribcell.cpp b/src/tallies/filter_distribcell.cpp
index c754dbd44ab..821e843daeb 100644
--- a/src/tallies/filter_distribcell.cpp
+++ b/src/tallies/filter_distribcell.cpp
@@ -1,5 +1,7 @@
 #include "openmc/tallies/filter_distribcell.h"
 
+#include <cassert>
+
 #include <fmt/core.h>
 
 #include "openmc/cell.h"
@@ -29,8 +31,8 @@ void DistribcellFilter::from_xml(pugi::xml_node node)
 
 void DistribcellFilter::set_cell(int32_t cell)
 {
-  Expects(cell >= 0);
-  Expects(cell < model::cells.size());
+  assert(cell >= 0);
+  assert(cell < model::cells.size());
   cell_ = cell;
   n_bins_ = model::cells[cell]->n_instances_;
 }
diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp
index 4767dd175f0..0b954cce3a3 100644
--- a/src/tallies/filter_energy.cpp
+++ b/src/tallies/filter_energy.cpp
@@ -21,14 +21,14 @@ void EnergyFilter::from_xml(pugi::xml_node node)
   this->set_bins(bins);
 }
 
-void EnergyFilter::set_bins(gsl::span<const double> bins)
+void EnergyFilter::set_bins(span<const double> bins)
 {
   // Clear existing bins
   bins_.clear();
   bins_.reserve(bins.size());
 
   // Copy bins, ensuring they are valid
-  for (gsl::index i = 0; i < bins.size(); ++i) {
+  for (int64_t i = 0; i < bins.size(); ++i) {
     if (i > 0 && bins[i] <= bins[i - 1]) {
       throw std::runtime_error {
         "Energy bins must be monotonically increasing."};
@@ -46,7 +46,7 @@ void EnergyFilter::set_bins(gsl::span<const double> bins)
   if (!settings::run_CE) {
     if (n_bins_ == data::mg.num_energy_groups_) {
       matches_transport_groups_ = true;
-      for (gsl::index i = 0; i < n_bins_ + 1; ++i) {
+      for (int64_t i = 0; i < n_bins_ + 1; ++i) {
         if (data::mg.rev_energy_bins_[i] != bins_[i]) {
           matches_transport_groups_ = false;
           break;
diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp
index 93ae24b2ada..fc4ba0ff91e 100644
--- a/src/tallies/filter_energyfunc.cpp
+++ b/src/tallies/filter_energyfunc.cpp
@@ -36,7 +36,7 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node)
 }
 
 void EnergyFunctionFilter::set_data(
-  gsl::span<const double> energy, gsl::span<const double> y)
+  span<const double> energy, span<const double> y)
 {
   // Check for consistent sizes with new data
   if (energy.size() != y.size()) {
@@ -48,7 +48,7 @@ void EnergyFunctionFilter::set_data(
   y_.reserve(y.size());
 
   // Copy over energy values, ensuring they are valid
-  for (gsl::index i = 0; i < energy.size(); ++i) {
+  for (int64_t i = 0; i < energy.size(); ++i) {
     if (i > 0 && energy[i] <= energy[i - 1]) {
       throw std::runtime_error {
         "Energy bins must be monotonically increasing."};
diff --git a/src/tallies/filter_material.cpp b/src/tallies/filter_material.cpp
index 6669ab2fa39..215c9af72ff 100644
--- a/src/tallies/filter_material.cpp
+++ b/src/tallies/filter_material.cpp
@@ -1,5 +1,7 @@
 #include "openmc/tallies/filter_material.h"
 
+#include <cassert>
+
 #include <fmt/core.h>
 
 #include "openmc/capi.h"
@@ -24,7 +26,7 @@ void MaterialFilter::from_xml(pugi::xml_node node)
   this->set_materials(mats);
 }
 
-void MaterialFilter::set_materials(gsl::span<const int32_t> materials)
+void MaterialFilter::set_materials(span<const int32_t> materials)
 {
   // Clear existing materials
   materials_.clear();
@@ -33,8 +35,8 @@ void MaterialFilter::set_materials(gsl::span<const int32_t> materials)
 
   // Update materials and mapping
   for (auto& index : materials) {
-    Expects(index >= 0);
-    Expects(index < model::materials.size());
+    assert(index >= 0);
+    assert(index < model::materials.size());
     materials_.push_back(index);
     map_[index] = materials_.size() - 1;
   }
diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp
index 03f7da97847..4edfbec4b9b 100644
--- a/src/tallies/filter_mesh.cpp
+++ b/src/tallies/filter_mesh.cpp
@@ -1,7 +1,6 @@
 #include "openmc/tallies/filter_mesh.h"
 
 #include <fmt/core.h>
-#include <gsl/gsl-lite.hpp>
 
 #include "openmc/capi.h"
 #include "openmc/constants.h"
diff --git a/src/tallies/filter_mu.cpp b/src/tallies/filter_mu.cpp
index 95bb3b21082..63915a53347 100644
--- a/src/tallies/filter_mu.cpp
+++ b/src/tallies/filter_mu.cpp
@@ -31,14 +31,14 @@ void MuFilter::from_xml(pugi::xml_node node)
   this->set_bins(bins);
 }
 
-void MuFilter::set_bins(gsl::span<double> bins)
+void MuFilter::set_bins(span<double> bins)
 {
   // Clear existing bins
   bins_.clear();
   bins_.reserve(bins.size());
 
   // Copy bins, ensuring they are valid
-  for (gsl::index i = 0; i < bins.size(); ++i) {
+  for (int64_t i = 0; i < bins.size(); ++i) {
     if (i > 0 && bins[i] <= bins[i - 1]) {
       throw std::runtime_error {"Mu bins must be monotonically increasing."};
     }
diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp
index 56d93a47f40..eef1d1e63c5 100644
--- a/src/tallies/filter_particle.cpp
+++ b/src/tallies/filter_particle.cpp
@@ -18,7 +18,7 @@ void ParticleFilter::from_xml(pugi::xml_node node)
   this->set_particles(types);
 }
 
-void ParticleFilter::set_particles(gsl::span<ParticleType> particles)
+void ParticleFilter::set_particles(span<ParticleType> particles)
 {
   // Clear existing particles
   particles_.clear();
diff --git a/src/tallies/filter_polar.cpp b/src/tallies/filter_polar.cpp
index d132ccf4236..29be6a437e3 100644
--- a/src/tallies/filter_polar.cpp
+++ b/src/tallies/filter_polar.cpp
@@ -32,14 +32,14 @@ void PolarFilter::from_xml(pugi::xml_node node)
   this->set_bins(bins);
 }
 
-void PolarFilter::set_bins(gsl::span<double> bins)
+void PolarFilter::set_bins(span<double> bins)
 {
   // Clear existing bins
   bins_.clear();
   bins_.reserve(bins.size());
 
   // Copy bins, ensuring they are valid
-  for (gsl::index i = 0; i < bins.size(); ++i) {
+  for (int64_t i = 0; i < bins.size(); ++i) {
     if (i > 0 && bins[i] <= bins[i - 1]) {
       throw std::runtime_error {"Polar bins must be monotonically increasing."};
     }
diff --git a/src/tallies/filter_sph_harm.cpp b/src/tallies/filter_sph_harm.cpp
index 359df379b2f..7441f32fb88 100644
--- a/src/tallies/filter_sph_harm.cpp
+++ b/src/tallies/filter_sph_harm.cpp
@@ -1,9 +1,9 @@
 #include "openmc/tallies/filter_sph_harm.h"
 
+#include <cassert>
 #include <utility> // For pair
 
 #include <fmt/core.h>
-#include <gsl/gsl-lite.hpp>
 
 #include "openmc/capi.h"
 #include "openmc/error.h"
@@ -30,7 +30,7 @@ void SphericalHarmonicsFilter::set_order(int order)
   n_bins_ = (order_ + 1) * (order_ + 1);
 }
 
-void SphericalHarmonicsFilter::set_cosine(gsl::cstring_span cosine)
+void SphericalHarmonicsFilter::set_cosine(const std::string& cosine)
 {
   if (cosine == "scatter") {
     cosine_ = SphericalHarmonicsCosine::scatter;
@@ -39,7 +39,7 @@ void SphericalHarmonicsFilter::set_cosine(gsl::cstring_span cosine)
   } else {
     throw std::invalid_argument {fmt::format("Unrecognized cosine type, \"{}\" "
                                              "in spherical harmonics filter",
-      gsl::to_string(cosine))};
+      cosine)};
   }
 }
 
@@ -87,7 +87,7 @@ void SphericalHarmonicsFilter::to_statepoint(hid_t filter_group) const
 
 std::string SphericalHarmonicsFilter::text_label(int bin) const
 {
-  Expects(bin >= 0 && bin < n_bins_);
+  assert(bin >= 0 && bin < n_bins_);
   for (int n = 0; n < order_ + 1; n++) {
     if (bin < (n + 1) * (n + 1)) {
       int m = (bin - n * n) - n;
diff --git a/src/tallies/filter_surface.cpp b/src/tallies/filter_surface.cpp
index a84d1772822..839a665c96a 100644
--- a/src/tallies/filter_surface.cpp
+++ b/src/tallies/filter_surface.cpp
@@ -1,5 +1,7 @@
 #include "openmc/tallies/filter_surface.h"
 
+#include <cassert>
+
 #include <fmt/core.h>
 
 #include "openmc/error.h"
@@ -26,7 +28,7 @@ void SurfaceFilter::from_xml(pugi::xml_node node)
   this->set_surfaces(surfaces);
 }
 
-void SurfaceFilter::set_surfaces(gsl::span<int32_t> surfaces)
+void SurfaceFilter::set_surfaces(span<int32_t> surfaces)
 {
   // Clear existing surfaces
   surfaces_.clear();
@@ -35,8 +37,8 @@ void SurfaceFilter::set_surfaces(gsl::span<int32_t> surfaces)
 
   // Update surfaces and mapping
   for (auto& index : surfaces) {
-    Expects(index >= 0);
-    Expects(index < model::surfaces.size());
+    assert(index >= 0);
+    assert(index < model::surfaces.size());
     surfaces_.push_back(index);
     map_[index] = surfaces_.size() - 1;
   }
diff --git a/src/tallies/filter_time.cpp b/src/tallies/filter_time.cpp
index 0378aa28de3..948d1347a0f 100644
--- a/src/tallies/filter_time.cpp
+++ b/src/tallies/filter_time.cpp
@@ -21,7 +21,7 @@ void TimeFilter::from_xml(pugi::xml_node node)
   this->set_bins(bins);
 }
 
-void TimeFilter::set_bins(gsl::span<const double> bins)
+void TimeFilter::set_bins(span<const double> bins)
 {
   // Clear existing bins
   bins_.clear();
diff --git a/src/tallies/filter_universe.cpp b/src/tallies/filter_universe.cpp
index 96ad0212d50..48114cfad59 100644
--- a/src/tallies/filter_universe.cpp
+++ b/src/tallies/filter_universe.cpp
@@ -1,5 +1,7 @@
 #include "openmc/tallies/filter_universe.h"
 
+#include <cassert>
+
 #include <fmt/core.h>
 
 #include "openmc/cell.h"
@@ -24,7 +26,7 @@ void UniverseFilter::from_xml(pugi::xml_node node)
   this->set_universes(universes);
 }
 
-void UniverseFilter::set_universes(gsl::span<int32_t> universes)
+void UniverseFilter::set_universes(span<int32_t> universes)
 {
   // Clear existing universes
   universes_.clear();
@@ -33,8 +35,8 @@ void UniverseFilter::set_universes(gsl::span<int32_t> universes)
 
   // Update universes and mapping
   for (auto& index : universes) {
-    Expects(index >= 0);
-    Expects(index < model::universes.size());
+    assert(index >= 0);
+    assert(index < model::universes.size());
     universes_.push_back(index);
     map_[index] = universes_.size() - 1;
   }
diff --git a/src/tallies/filter_zernike.cpp b/src/tallies/filter_zernike.cpp
index eb4c8bdfd00..af5b595aabb 100644
--- a/src/tallies/filter_zernike.cpp
+++ b/src/tallies/filter_zernike.cpp
@@ -1,11 +1,11 @@
 #include "openmc/tallies/filter_zernike.h"
 
+#include <cassert>
 #include <cmath>
 #include <sstream>
 #include <utility> // For pair
 
 #include <fmt/core.h>
-#include <gsl/gsl-lite.hpp>
 
 #include "openmc/capi.h"
 #include "openmc/error.h"
@@ -57,7 +57,7 @@ void ZernikeFilter::to_statepoint(hid_t filter_group) const
 
 std::string ZernikeFilter::text_label(int bin) const
 {
-  Expects(bin >= 0 && bin < n_bins_);
+  assert(bin >= 0 && bin < n_bins_);
   for (int n = 0; n < order_ + 1; n++) {
     int last = (n + 1) * (n + 2) / 2;
     if (bin < last) {
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 4f33abf6beb..96d684f71a8 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -39,7 +39,8 @@
 #include <fmt/core.h>
 
 #include <algorithm> // for max
-#include <cstddef>   // for size_t
+#include <cassert>
+#include <cstddef> // for size_t
 #include <string>
 
 namespace openmc {
@@ -140,7 +141,7 @@ Tally::Tally(pugi::xml_node node)
   // Check for the presence of certain filter types
   bool has_energyout = energyout_filter_ >= 0;
   int particle_filter_index = C_NONE;
-  for (gsl::index j = 0; j < filters_.size(); ++j) {
+  for (int64_t j = 0; j < filters_.size(); ++j) {
     int i_filter = filters_[j];
     const auto& f = model::tally_filters[i_filter].get();
 
@@ -345,7 +346,7 @@ Tally* Tally::create(int32_t id)
 
 void Tally::set_id(int32_t id)
 {
-  Expects(id >= 0 || id == C_NONE);
+  assert(id >= 0 || id == C_NONE);
 
   // Clear entry in tally map if an ID was already assigned before
   if (id_ != C_NONE) {
@@ -401,7 +402,7 @@ bool Tally::has_filter(FilterType filter_type) const
   return false;
 }
 
-void Tally::set_filters(gsl::span<Filter*> filters)
+void Tally::set_filters(span<Filter*> filters)
 {
   // Clear old data.
   filters_.clear();
@@ -1370,7 +1371,7 @@ extern "C" int openmc_tally_set_filters(
   try {
     // Convert indices to filter pointers
     vector<Filter*> filters;
-    for (gsl::index i = 0; i < n; ++i) {
+    for (int64_t i = 0; i < n; ++i) {
       int32_t i_filt = indices[i];
       filters.push_back(model::tally_filters.at(i_filt).get());
     }
diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp
index 68f7550ae5a..beab9e50d1f 100644
--- a/src/weight_windows.cpp
+++ b/src/weight_windows.cpp
@@ -1,6 +1,7 @@
 #include "openmc/weight_windows.h"
 
 #include <algorithm>
+#include <cassert>
 #include <cmath>
 #include <set>
 #include <string>
@@ -31,7 +32,6 @@
 #include "openmc/xml_interface.h"
 
 #include <fmt/core.h>
-#include <gsl/gsl-lite.hpp>
 
 namespace openmc {
 
@@ -288,7 +288,7 @@ void WeightWindows::allocate_ww_bounds()
 
 void WeightWindows::set_id(int32_t id)
 {
-  Expects(id >= 0 || id == C_NONE);
+  assert(id >= 0 || id == C_NONE);
 
   // Clear entry in mesh map in case one was already assigned
   if (id_ != C_NONE) {
@@ -316,7 +316,7 @@ void WeightWindows::set_id(int32_t id)
   variance_reduction::ww_map[id] = index_;
 }
 
-void WeightWindows::set_energy_bounds(gsl::span<const double> bounds)
+void WeightWindows::set_energy_bounds(span<const double> bounds)
 {
   energy_bounds_.clear();
   energy_bounds_.insert(energy_bounds_.begin(), bounds.begin(), bounds.end());
@@ -451,7 +451,7 @@ void WeightWindows::set_bounds(
 }
 
 void WeightWindows::set_bounds(
-  gsl::span<const double> lower_bounds, gsl::span<const double> upper_bounds)
+  span<const double> lower_bounds, span<const double> upper_bounds)
 {
   check_bounds(lower_bounds, upper_bounds);
   auto shape = this->bounds_size();
@@ -465,8 +465,7 @@ void WeightWindows::set_bounds(
     xt::adapt(upper_bounds.data(), upper_ww_.shape());
 }
 
-void WeightWindows::set_bounds(
-  gsl::span<const double> lower_bounds, double ratio)
+void WeightWindows::set_bounds(span<const double> lower_bounds, double ratio)
 {
   this->check_bounds(lower_bounds);
 
diff --git a/vendor/gsl-lite b/vendor/gsl-lite
deleted file mode 160000
index 913e86d49c6..00000000000
--- a/vendor/gsl-lite
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit 913e86d49c6a1acca980f4e325378f9dc393493a