diff --git a/src/axom/core/CMakeLists.txt b/src/axom/core/CMakeLists.txt
index fc4be76825..39627d95ac 100644
--- a/src/axom/core/CMakeLists.txt
+++ b/src/axom/core/CMakeLists.txt
@@ -60,6 +60,7 @@ set(core_headers
     Macros.hpp
     Map.hpp
     FlatMap.hpp
+    NumericLimits.hpp
     Path.hpp
     RangeAdapter.hpp
     StackArray.hpp
diff --git a/src/axom/core/NumericLimits.hpp b/src/axom/core/NumericLimits.hpp
new file mode 100644
index 0000000000..8e5642b8fa
--- /dev/null
+++ b/src/axom/core/NumericLimits.hpp
@@ -0,0 +1,41 @@
+// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and
+// other Axom Project Developers. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+/*!
+ *
+ * \file NumericLimits.hpp
+ *
+ * \brief Header file containing portability layer for std::numeric_limits
+ *        capabilities
+ *
+ */
+
+#ifndef AXOM_NUMERICLIMITS_HPP_
+#define AXOM_NUMERICLIMITS_HPP_
+
+#include "axom/config.hpp"  // for compile-time definitions
+
+#include <limits>
+
+#if defined(AXOM_USE_CUDA)
+  #include <cuda/std/limits>
+#endif
+
+namespace axom
+{
+#if defined(AXOM_USE_CUDA)
+// Note: cuda::std types work in host and device code as long as Axom is
+//       configured with CUDA enabled. No need to rely on two different
+//       header files in that case.
+template <typename T>
+using numeric_limits = cuda::std::numeric_limits<T>;
+#else
+template <typename T>
+using numeric_limits = std::numeric_limits<T>;
+#endif
+
+}  // namespace axom
+
+#endif  // AXOM_NUMERICLIMITS_HPP_
diff --git a/src/axom/core/tests/CMakeLists.txt b/src/axom/core/tests/CMakeLists.txt
index bf77f9951e..d1ebe2810d 100644
--- a/src/axom/core/tests/CMakeLists.txt
+++ b/src/axom/core/tests/CMakeLists.txt
@@ -26,6 +26,7 @@ set(core_serial_tests
     core_map.hpp
     core_flatmap.hpp
     core_memory_management.hpp
+    core_numeric_limits.hpp
     core_Path.hpp
     core_stack_array.hpp
     core_static_array.hpp
diff --git a/src/axom/core/tests/core_bit_utilities.hpp b/src/axom/core/tests/core_bit_utilities.hpp
index af3cc03934..d8ec6eb5e4 100644
--- a/src/axom/core/tests/core_bit_utilities.hpp
+++ b/src/axom/core/tests/core_bit_utilities.hpp
@@ -7,6 +7,7 @@
 
 #include "axom/config.hpp"
 #include "axom/core/Types.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/core/utilities/Utilities.hpp"
 #include "axom/core/utilities/BitUtilities.hpp"
 
@@ -22,7 +23,7 @@ T random_int()
 {
   static_assert(std::is_integral<T>::value, "T must be an integral type");
 
-  constexpr T max_int = std::numeric_limits<T>::max();
+  constexpr T max_int = axom::numeric_limits<T>::max();
   constexpr double max_d = static_cast<double>(max_int);
 
   const auto val = axom::utilities::random_real(0., max_d);
diff --git a/src/axom/core/tests/core_numeric_limits.hpp b/src/axom/core/tests/core_numeric_limits.hpp
new file mode 100644
index 0000000000..18727579cc
--- /dev/null
+++ b/src/axom/core/tests/core_numeric_limits.hpp
@@ -0,0 +1,170 @@
+// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and
+// other Axom Project Developers. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#include "axom/config.hpp"  // for compile time definitions
+
+#include "axom/core/NumericLimits.hpp"
+
+// for gtest macros
+#include "gtest/gtest.h"
+
+//------------------------------------------------------------------------------
+//  UNIT TESTS
+//------------------------------------------------------------------------------
+
+//------------------------------------------------------------------------------
+TEST(core_NumericLimits, check_CPU)
+{
+  //
+  // Tests to compare axom::numeric_limits to std::numeric_limits
+  // to ensure that Axom type aliasing is correct.
+  //
+  EXPECT_TRUE(axom::numeric_limits<int>::lowest() ==
+              std::numeric_limits<int>::lowest());
+  EXPECT_TRUE(axom::numeric_limits<int>::min() == std::numeric_limits<int>::min());
+  EXPECT_TRUE(axom::numeric_limits<int>::max() == std::numeric_limits<int>::max());
+  EXPECT_TRUE(axom::numeric_limits<int>::is_signed ==
+              std::numeric_limits<int>::is_signed);
+
+  EXPECT_TRUE(axom::numeric_limits<float>::lowest() ==
+              std::numeric_limits<float>::lowest());
+  EXPECT_TRUE(axom::numeric_limits<float>::min() ==
+              std::numeric_limits<float>::min());
+  EXPECT_TRUE(axom::numeric_limits<float>::max() ==
+              std::numeric_limits<float>::max());
+
+  EXPECT_TRUE(axom::numeric_limits<double>::lowest() ==
+              std::numeric_limits<double>::lowest());
+  EXPECT_TRUE(axom::numeric_limits<double>::min() ==
+              std::numeric_limits<double>::min());
+  EXPECT_TRUE(axom::numeric_limits<double>::max() ==
+              std::numeric_limits<double>::max());
+}
+
+//------------------------------------------------------------------------------
+#if defined(AXOM_USE_CUDA)
+//
+// Tests to ensure axom::numeric_limits type alias does the correct thing
+// in host and CUDA device code.
+//
+
+//
+// Simple device kernel
+//
+__global__ void cuda_kernel(int* a, size_t* b, float* c, double* d)
+{
+  a[0] = axom::numeric_limits<int>::min();
+  b[0] = axom::numeric_limits<size_t>::max();
+  c[0] = axom::numeric_limits<float>::lowest();
+  d[0] = axom::numeric_limits<double>::max();
+}
+
+TEST(core_NumericLimits, check_CUDA)
+{
+  //
+  // Device memory allocation and initialiation for a few different types.
+  //
+  int* a;
+  (void)cudaMalloc(&a, sizeof(int));
+  (void)cudaMemset(a, 0, sizeof(int));
+
+  size_t* b;
+  (void)cudaMalloc(&b, sizeof(size_t));
+  (void)cudaMemset(b, 0, sizeof(size_t));
+
+  float* c;
+  (void)cudaMalloc(&c, sizeof(float));
+  (void)cudaMemset(c, 0, sizeof(float));
+
+  double* d;
+  (void)cudaMalloc(&d, sizeof(double));
+  (void)cudaMemset(d, 0, sizeof(double));
+
+  //
+  // Set values in device code.
+  //
+  cuda_kernel<<<1, 1>>>(a, b, c, d);
+
+  //
+  // Copy device values back to host and compare with expectations....
+  //
+  int ha;
+  size_t hb;
+  float hc;
+  double hd;
+  (void)cudaMemcpy(&ha, a, sizeof(int), cudaMemcpyDeviceToHost);
+  (void)cudaMemcpy(&hb, b, sizeof(size_t), cudaMemcpyDeviceToHost);
+  (void)cudaMemcpy(&hc, c, sizeof(float), cudaMemcpyDeviceToHost);
+  (void)cudaMemcpy(&hd, d, sizeof(double), cudaMemcpyDeviceToHost);
+
+  EXPECT_TRUE(ha == axom::numeric_limits<int>::min());
+  EXPECT_TRUE(hb == axom::numeric_limits<size_t>::max());
+  EXPECT_TRUE(hc == axom::numeric_limits<float>::lowest());
+  EXPECT_TRUE(hd == axom::numeric_limits<double>::max());
+}
+#endif
+
+//------------------------------------------------------------------------------
+#if defined(AXOM_USE_HIP)
+//
+// Tests to ensure axom::numeric_limits type alias does the correct thing
+// in host and CUDA device code.
+//
+
+//
+// Simple device kernel
+//
+__global__ void hip_kernel(int* a, size_t* b, float* c, double* d)
+{
+  a[0] = axom::numeric_limits<int>::min();
+  b[0] = axom::numeric_limits<size_t>::max();
+  c[0] = axom::numeric_limits<float>::lowest();
+  d[0] = axom::numeric_limits<double>::max();
+}
+
+TEST(core_NumericLimits, check_HIP)
+{
+  //
+  // Device memory allocation and initialiation for a few different types.
+  //
+  int* a;
+  (void)hipMalloc(&a, sizeof(int));
+  (void)hipMemset(a, 0, sizeof(int));
+
+  size_t* b;
+  (void)hipMalloc(&b, sizeof(size_t));
+  (void)hipMemset(b, 0, sizeof(size_t));
+
+  float* c;
+  (void)hipMalloc(&c, sizeof(float));
+  (void)hipMemset(c, 0, sizeof(float));
+
+  double* d;
+  (void)hipMalloc(&d, sizeof(double));
+  (void)hipMemset(d, 0, sizeof(double));
+
+  //
+  // Set values in device code.
+  //
+  hip_kernel<<<1, 1>>>(a, b, c, d);
+
+  //
+  // Copy device values back to host and compare with expectations....
+  //
+  int ha;
+  size_t hb;
+  float hc;
+  double hd;
+  (void)hipMemcpy(&ha, a, sizeof(int), hipMemcpyDeviceToHost);
+  (void)hipMemcpy(&hb, b, sizeof(size_t), hipMemcpyDeviceToHost);
+  (void)hipMemcpy(&hc, c, sizeof(float), hipMemcpyDeviceToHost);
+  (void)hipMemcpy(&hd, d, sizeof(double), hipMemcpyDeviceToHost);
+
+  EXPECT_TRUE(ha == axom::numeric_limits<int>::min());
+  EXPECT_TRUE(hb == axom::numeric_limits<size_t>::max());
+  EXPECT_TRUE(hc == axom::numeric_limits<float>::lowest());
+  EXPECT_TRUE(hd == axom::numeric_limits<double>::max());
+}
+#endif
diff --git a/src/axom/core/tests/core_serial_main.cpp b/src/axom/core/tests/core_serial_main.cpp
index 060318d31e..659a607861 100644
--- a/src/axom/core/tests/core_serial_main.cpp
+++ b/src/axom/core/tests/core_serial_main.cpp
@@ -18,6 +18,7 @@
 #include "core_map.hpp"
 #include "core_flatmap.hpp"
 #include "core_memory_management.hpp"
+#include "core_numeric_limits.hpp"
 #include "core_Path.hpp"
 #include "core_stack_array.hpp"
 #include "core_static_array.hpp"
diff --git a/src/axom/core/tests/core_types.hpp b/src/axom/core/tests/core_types.hpp
index 6a5cfd4ce1..d60b2d4486 100644
--- a/src/axom/core/tests/core_types.hpp
+++ b/src/axom/core/tests/core_types.hpp
@@ -7,12 +7,12 @@
 #include "axom/config.hpp"
 #include "axom/core/Types.hpp"
 #include "axom/core/Macros.hpp"
+#include "axom/core/NumericLimits.hpp"
 
 // gtest includes
 #include "gtest/gtest.h"
 
 // C/C++ includes
-#include <limits>       // for std::numeric_limits
 #include <type_traits>  // for std::is_same, std::is_integral, etc.
 
 #ifndef AXOM_USE_MPI
@@ -62,7 +62,7 @@ void check_real_type(std::size_t expected_num_bytes,
                      MPI_Datatype expected_mpi_type)
 {
   EXPECT_TRUE(std::is_floating_point<RealType>::value);
-  EXPECT_TRUE(std::numeric_limits<RealType>::is_signed);
+  EXPECT_TRUE(axom::numeric_limits<RealType>::is_signed);
   EXPECT_EQ(sizeof(RealType), expected_num_bytes);
 
   check_mpi_type<RealType>(expected_num_bytes, expected_mpi_type);
@@ -75,9 +75,9 @@ void check_integral_type(std::size_t expected_num_bytes,
                          int expected_num_digits,
                          MPI_Datatype expected_mpi_type)
 {
-  EXPECT_TRUE(std::numeric_limits<IntegralType>::is_integer);
-  EXPECT_EQ(std::numeric_limits<IntegralType>::is_signed, is_signed);
-  EXPECT_EQ(std::numeric_limits<IntegralType>::digits, expected_num_digits);
+  EXPECT_TRUE(axom::numeric_limits<IntegralType>::is_integer);
+  EXPECT_EQ(axom::numeric_limits<IntegralType>::is_signed, is_signed);
+  EXPECT_EQ(axom::numeric_limits<IntegralType>::digits, expected_num_digits);
   EXPECT_EQ(sizeof(IntegralType), expected_num_bytes);
 
   check_mpi_type<IntegralType>(expected_num_bytes, expected_mpi_type);
diff --git a/src/axom/core/tests/numerics_floating_point_limits.hpp b/src/axom/core/tests/numerics_floating_point_limits.hpp
index a9e4bab880..833c570086 100644
--- a/src/axom/core/tests/numerics_floating_point_limits.hpp
+++ b/src/axom/core/tests/numerics_floating_point_limits.hpp
@@ -47,5 +47,7 @@ TEST(numerics_floating_point_limits, consistency_with_standard_numeric_limits)
 {
   check_type_limits<float>("float");
   check_type_limits<double>("double");
+#if !defined(AXOM_DEVICE_CODE)
   check_type_limits<long double>("long double");
+#endif
 }
diff --git a/src/axom/inlet/examples/documentation_generation.cpp b/src/axom/inlet/examples/documentation_generation.cpp
index c9748b4cf6..ed74434d0d 100644
--- a/src/axom/inlet/examples/documentation_generation.cpp
+++ b/src/axom/inlet/examples/documentation_generation.cpp
@@ -6,11 +6,11 @@
 // usage : ./inlet_documentation_generation_example --enableDocs --fil lua_file.lua
 
 #include "axom/inlet.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/slic/core/SimpleLogger.hpp"
 
 #include "axom/CLI11.hpp"
 #include <iostream>
-#include <limits>
 
 using axom::inlet::Inlet;
 using axom::inlet::LuaReader;
@@ -68,17 +68,17 @@ void defineSchema(Inlet& inlet)
   filename_field.required();
 
   inlet.addInt("thermal_solver/mesh/serial", "number of serial refinements")
-    .range(0, std::numeric_limits<int>::max())
+    .range(0, axom::numeric_limits<int>::max())
     .defaultValue(1);
 
   // The description for thermal_solver/mesh/parallel is left unspecified
   inlet.addInt("thermal_solver/mesh/parallel")
-    .range(1, std::numeric_limits<int>::max())
+    .range(1, axom::numeric_limits<int>::max())
     .defaultValue(1);
 
   inlet.addInt("thermal_solver/order", "polynomial order")
     .required()
-    .range(1, std::numeric_limits<int>::max());
+    .range(1, axom::numeric_limits<int>::max());
 
   auto& timestep_field =
     inlet.addString("thermal_solver/timestepper", "thermal solver timestepper");
@@ -109,13 +109,13 @@ void defineSchema(Inlet& inlet)
     solver_schema.addDouble("rel_tol", "solver relative tolerance");
   rel_tol_field.required(false);
   rel_tol_field.defaultValue(1.e-6);
-  rel_tol_field.range(0.0, std::numeric_limits<double>::max());
+  rel_tol_field.range(0.0, axom::numeric_limits<double>::max());
 
   auto& abs_tol_field =
     solver_schema.addDouble("abs_tol", "solver absolute tolerance");
   abs_tol_field.required(true);
   abs_tol_field.defaultValue(1.e-12);
-  abs_tol_field.range(0.0, std::numeric_limits<double>::max());
+  abs_tol_field.range(0.0, axom::numeric_limits<double>::max());
 
   auto& print_level_field =
     solver_schema.addInt("print_level", "solver print/debug level");
@@ -127,18 +127,18 @@ void defineSchema(Inlet& inlet)
     solver_schema.addInt("max_iter", "maximum iteration limit");
   max_iter_field.required(false);
   max_iter_field.defaultValue(100);
-  max_iter_field.range(1, std::numeric_limits<int>::max());
+  max_iter_field.range(1, axom::numeric_limits<int>::max());
 
   auto& dt_field = solver_schema.addDouble("dt", "time step");
   dt_field.required(true);
   dt_field.defaultValue(1);
-  dt_field.range(0.0, std::numeric_limits<double>::max());
+  dt_field.range(0.0, axom::numeric_limits<double>::max());
 
   auto& steps_field =
     solver_schema.addInt("steps", "number of steps/cycles to take");
   steps_field.required(true);
   steps_field.defaultValue(1);
-  steps_field.range(1, std::numeric_limits<int>::max());
+  steps_field.range(1, axom::numeric_limits<int>::max());
 }
 
 // Checking the contents of the passed inlet
diff --git a/src/axom/mint/mesh/StructuredMesh.cpp b/src/axom/mint/mesh/StructuredMesh.cpp
index 067028b821..2550a92876 100644
--- a/src/axom/mint/mesh/StructuredMesh.cpp
+++ b/src/axom/mint/mesh/StructuredMesh.cpp
@@ -7,7 +7,7 @@
 #include "axom/mint/mesh/blueprint.hpp"
 
 #include <cstring> /* for memcpy() */
-#include <limits>  /* for std::numeric_limits< IndexType >::max */
+#include "axom/core/NumericLimits.hpp"
 
 namespace axom
 {
@@ -169,9 +169,9 @@ void StructuredMesh::structuredInit()
 
   /* Initialize the node meta data. */
   m_node_jp = (m_ndims > 1) ? getNodeResolution(0)
-                            : std::numeric_limits<IndexType>::max();
+                            : axom::numeric_limits<IndexType>::max();
   m_node_kp = (m_ndims > 2) ? m_node_jp * getNodeResolution(1)
-                            : std::numeric_limits<IndexType>::max();
+                            : axom::numeric_limits<IndexType>::max();
 
   /* Initialize the cell meta data */
   for(int dim = 0; dim < m_ndims; ++dim)
@@ -180,9 +180,9 @@ void StructuredMesh::structuredInit()
   }
 
   m_cell_jp = (m_ndims > 1) ? getCellResolution(0)
-                            : std::numeric_limits<IndexType>::max();
+                            : axom::numeric_limits<IndexType>::max();
   m_cell_kp = (m_ndims > 2) ? m_cell_jp * getCellResolution(1)
-                            : std::numeric_limits<IndexType>::max();
+                            : axom::numeric_limits<IndexType>::max();
 
   /* Build the cell to node offsets. */
   m_cell_node_offsets[0] = 0;
diff --git a/src/axom/mint/tests/mint_fem_shape_functions.cpp b/src/axom/mint/tests/mint_fem_shape_functions.cpp
index ab1ea065b3..05ed8bdb65 100644
--- a/src/axom/mint/tests/mint_fem_shape_functions.cpp
+++ b/src/axom/mint/tests/mint_fem_shape_functions.cpp
@@ -5,6 +5,7 @@
 
 #include "gtest/gtest.h"
 
+#include "axom/core/NumericLimits.hpp"
 #include "axom/mint/fem/shape_functions/Lagrange.hpp"
 #include "axom/mint/fem/shape_functions/ShapeFunction.hpp"
 #include "axom/mint/fem/FEBasis.hpp"
@@ -13,9 +14,6 @@
 
 #include "axom/slic.hpp"
 
-// C/C++ includes
-#include <limits>
-
 using namespace axom;
 using mint::Lagrange;
 using mint::ShapeFunction;
@@ -32,7 +30,7 @@ namespace
  * \tparam CELLTYPE the corresponding cell type, e.g., MINT_QUAD
  */
 template <int BasisType, mint::CellType CELLTYPE>
-void reference_element(double TOL = std::numeric_limits<double>::epsilon())
+void reference_element(double TOL = axom::numeric_limits<double>::epsilon())
 {
   using FEMType = typename mint::FEBasis<BasisType, CELLTYPE>;
   using ShapeFunctionType = typename FEMType::ShapeFunctionType;
diff --git a/src/axom/primal/geometry/BoundingBox.hpp b/src/axom/primal/geometry/BoundingBox.hpp
index c767b69fe8..5f8e3a034b 100644
--- a/src/axom/primal/geometry/BoundingBox.hpp
+++ b/src/axom/primal/geometry/BoundingBox.hpp
@@ -6,12 +6,11 @@
 #ifndef AXOM_PRIMAL_BOUNDINGBOX_HPP_
 #define AXOM_PRIMAL_BOUNDINGBOX_HPP_
 
-#include <limits>
-
 #include "axom/config.hpp"
 
 #include "axom/core/Macros.hpp"
 #include "axom/core/numerics/floating_point_limits.hpp"
+#include "axom/core/NumericLimits.hpp"
 
 #include "axom/primal/geometry/Point.hpp"
 #include "axom/primal/geometry/Vector.hpp"
@@ -72,8 +71,8 @@ class BoundingBox
   using VectorType = Vector<T, NDIMS>;
   using BoxType = BoundingBox<T, NDIMS>;
 
-  static constexpr T InvalidMin = std::numeric_limits<T>::max();
-  static constexpr T InvalidMax = std::numeric_limits<T>::lowest();
+  static constexpr T InvalidMin = axom::numeric_limits<T>::max();
+  static constexpr T InvalidMax = axom::numeric_limits<T>::lowest();
 
 public:
   /*!
diff --git a/src/axom/primal/geometry/OrientedBoundingBox.hpp b/src/axom/primal/geometry/OrientedBoundingBox.hpp
index d45e1a0a48..13e0d31683 100644
--- a/src/axom/primal/geometry/OrientedBoundingBox.hpp
+++ b/src/axom/primal/geometry/OrientedBoundingBox.hpp
@@ -10,6 +10,7 @@
 
 #include "axom/config.hpp"
 #include "axom/core.hpp"
+#include "axom/core/NumericLimits.hpp"
 
 #include "axom/primal/geometry/NumericArray.hpp"
 #include "axom/primal/geometry/Point.hpp"
@@ -490,7 +491,7 @@ OrientedBoundingBox<T, NDIMS>::OrientedBoundingBox(const OrientedBoundingBox& ot
 template <typename T, int NDIMS>
 void OrientedBoundingBox<T, NDIMS>::clear()
 {
-  (this->m_u[0])[0] = std::numeric_limits<T>::lowest();
+  (this->m_u[0])[0] = axom::numeric_limits<T>::lowest();
 }
 
 //------------------------------------------------------------------------------
@@ -699,7 +700,7 @@ Point<T, NDIMS> OrientedBoundingBox<T, NDIMS>::furthestPoint(const PointType& pt
 template <typename T, int NDIMS>
 bool OrientedBoundingBox<T, NDIMS>::isValid() const
 {
-  return !((this->m_u[0])[0] == std::numeric_limits<T>::lowest());
+  return !((this->m_u[0])[0] == axom::numeric_limits<T>::lowest());
 }
 
 //------------------------------------------------------------------------------
diff --git a/src/axom/primal/tests/primal_boundingbox.cpp b/src/axom/primal/tests/primal_boundingbox.cpp
index 64f862f9ab..1a15212c0d 100644
--- a/src/axom/primal/tests/primal_boundingbox.cpp
+++ b/src/axom/primal/tests/primal_boundingbox.cpp
@@ -3,13 +3,12 @@
 //
 // SPDX-License-Identifier: (BSD-3-Clause)
 
-#include <limits>
-
 #include "gtest/gtest.h"
 #include "axom/slic.hpp"
 
 #include "axom/core/execution/execution_space.hpp"
 #include "axom/core/execution/for_all.hpp"
+#include "axom/core/NumericLimits.hpp"
 
 #include "axom/primal/geometry/Point.hpp"
 #include "axom/primal/geometry/BoundingBox.hpp"
@@ -612,28 +611,28 @@ TEST(primal_boundingBox, highest_lowest_values)
   // is doing the right thing in our CXX11 and pre-CXX11 compilers
 
   // Test double
-  double maxD = std::numeric_limits<double>::max();
+  double maxD = axom::numeric_limits<double>::max();
   double minD = -maxD;
-  EXPECT_EQ(maxD, std::numeric_limits<double>::max());
-  EXPECT_EQ(minD, std::numeric_limits<double>::lowest());
+  EXPECT_EQ(maxD, axom::numeric_limits<double>::max());
+  EXPECT_EQ(minD, axom::numeric_limits<double>::lowest());
 
   // Test float
-  double maxF = std::numeric_limits<float>::max();
+  double maxF = axom::numeric_limits<float>::max();
   double minF = -maxF;
-  EXPECT_EQ(maxF, std::numeric_limits<float>::max());
-  EXPECT_EQ(minF, std::numeric_limits<float>::lowest());
+  EXPECT_EQ(maxF, axom::numeric_limits<float>::max());
+  EXPECT_EQ(minF, axom::numeric_limits<float>::lowest());
 
   // Test int
-  int maxI = std::numeric_limits<int>::max();
-  int minI = std::numeric_limits<int>::min();
-  EXPECT_EQ(maxI, std::numeric_limits<int>::max());
-  EXPECT_EQ(minI, std::numeric_limits<int>::lowest());
+  int maxI = axom::numeric_limits<int>::max();
+  int minI = axom::numeric_limits<int>::min();
+  EXPECT_EQ(maxI, axom::numeric_limits<int>::max());
+  EXPECT_EQ(minI, axom::numeric_limits<int>::lowest());
 
   // Test uint
-  unsigned int maxU = std::numeric_limits<unsigned int>::max();
-  unsigned int minU = std::numeric_limits<unsigned int>::min();
-  EXPECT_EQ(maxU, std::numeric_limits<unsigned int>::max());
-  EXPECT_EQ(minU, std::numeric_limits<unsigned int>::lowest());
+  unsigned int maxU = axom::numeric_limits<unsigned int>::max();
+  unsigned int minU = axom::numeric_limits<unsigned int>::min();
+  EXPECT_EQ(maxU, axom::numeric_limits<unsigned int>::max());
+  EXPECT_EQ(minU, axom::numeric_limits<unsigned int>::lowest());
 
   // Testing that our default constructor for bounding boxes is properly setting the range.
 
diff --git a/src/axom/primal/tests/primal_squared_distance.cpp b/src/axom/primal/tests/primal_squared_distance.cpp
index d780694a49..c9a3b77b88 100644
--- a/src/axom/primal/tests/primal_squared_distance.cpp
+++ b/src/axom/primal/tests/primal_squared_distance.cpp
@@ -6,6 +6,7 @@
 #include "gtest/gtest.h"
 
 #include "axom/config.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/slic.hpp"
 #include "axom/primal.hpp"
 
@@ -201,7 +202,7 @@ TEST(primal_squared_distance, point_to_bbox)
           EXPECT_NEAR(sqsum, primal::squared_distance(pt, cube), EPS);
         }
 
-        EXPECT_EQ(std::numeric_limits<double>::max(),
+        EXPECT_EQ(axom::numeric_limits<double>::max(),
                   squared_distance(pt, empty));
       }
     }
@@ -233,9 +234,9 @@ TEST(primal_squared_distance, bbox_to_bbox)
 
   // check that squared distances for empty/invalid boxes is max double
   const QBBox empty;
-  EXPECT_EQ(std::numeric_limits<double>::max(), squared_distance(middle, empty));
-  EXPECT_EQ(std::numeric_limits<double>::max(), squared_distance(empty, middle));
-  EXPECT_EQ(std::numeric_limits<double>::max(), squared_distance(empty, empty));
+  EXPECT_EQ(axom::numeric_limits<double>::max(), squared_distance(middle, empty));
+  EXPECT_EQ(axom::numeric_limits<double>::max(), squared_distance(empty, middle));
+  EXPECT_EQ(axom::numeric_limits<double>::max(), squared_distance(empty, empty));
 }
 
 //----------------------------------------------------------------------
diff --git a/src/axom/quest/DistributedClosestPoint.cpp b/src/axom/quest/DistributedClosestPoint.cpp
index 9bf6063a19..77c53b977f 100644
--- a/src/axom/quest/DistributedClosestPoint.cpp
+++ b/src/axom/quest/DistributedClosestPoint.cpp
@@ -8,6 +8,7 @@
 #include "axom/core/memory_management.hpp"
 #include "axom/core/execution/execution_space.hpp"
 #include "axom/core/execution/runtime_policy.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/quest/DistributedClosestPoint.hpp"
 #include "axom/quest/detail/DistributedClosestPointImpl.hpp"
 
@@ -17,7 +18,6 @@
 #include "conduit_blueprint_mpi.hpp"
 #include "conduit_relay_mpi.hpp"
 
-#include <limits>
 #include <cstdlib>
 
 #ifndef AXOM_USE_MPI
@@ -33,7 +33,7 @@ DistributedClosestPoint::DistributedClosestPoint()
   : m_mpiComm(MPI_COMM_WORLD)
   , m_mpiCommIsPrivate(false)
   , m_allocatorID(axom::INVALID_ALLOCATOR_ID)
-  , m_sqDistanceThreshold(std::numeric_limits<double>::max())
+  , m_sqDistanceThreshold(axom::numeric_limits<double>::max())
 {
   setDefaultAllocatorID();
   setMpiCommunicator(MPI_COMM_WORLD);
diff --git a/src/axom/quest/InOutOctree.hpp b/src/axom/quest/InOutOctree.hpp
index 93a8862de9..ebd99fdf94 100644
--- a/src/axom/quest/InOutOctree.hpp
+++ b/src/axom/quest/InOutOctree.hpp
@@ -13,6 +13,7 @@
 #define AXOM_QUEST_INOUT_OCTREE__HPP_
 
 #include "axom/core.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/slic.hpp"
 #include "axom/slam.hpp"
 #include "axom/primal.hpp"
@@ -28,7 +29,6 @@
 
 #include <vector>
 #include <iterator>
-#include <limits>
 #include <sstream>
 #include <unordered_map>
 
@@ -1203,7 +1203,7 @@ typename std::enable_if<TDIM == 3, bool>::type InOutOctree<DIM>::withinGrayBlock
     /// intersection. Note: We have to check all triangles to ensure that
     /// there is not a closer triangle than tri along this direction.
     CellIndex tIdx = MeshWrapper<DIM>::NO_CELL;
-    double minRayParam = std::numeric_limits<double>::infinity();
+    double minRayParam = axom::numeric_limits<double>::infinity();
     SpaceRay ray(queryPt, SpaceVector(queryPt, triPt));
 
     QUEST_OCTREE_DEBUG_LOG_IF(
@@ -1320,8 +1320,8 @@ typename std::enable_if<TDIM == 2, bool>::type InOutOctree<DIM>::withinGrayBlock
     // Using a ray from query pt to point on this segment
     // Find closest intersection to surface within cell inside this bounding box
     CellIndex tIdx = MeshWrapper<DIM>::NO_CELL;
-    double minRayParam = std::numeric_limits<double>::infinity();
-    double minSegParam = std::numeric_limits<double>::infinity();
+    double minRayParam = axom::numeric_limits<double>::infinity();
+    double minSegParam = axom::numeric_limits<double>::infinity();
     SpaceRay ray(queryPt, SpaceVector(queryPt, segmentPt));
 
     QUEST_OCTREE_DEBUG_LOG_IF(
diff --git a/src/axom/quest/MeshViewUtil.hpp b/src/axom/quest/MeshViewUtil.hpp
index 7c3098c6e5..4e38efcdd4 100644
--- a/src/axom/quest/MeshViewUtil.hpp
+++ b/src/axom/quest/MeshViewUtil.hpp
@@ -12,6 +12,7 @@
 #ifdef AXOM_USE_CONDUIT
 
   #include "axom/core.hpp"
+  #include "axom/core/NumericLimits.hpp"
   #include "axom/fmt.hpp"
   #include "conduit_blueprint.hpp"
   #include "conduit_blueprint_mcarray.hpp"
@@ -21,7 +22,6 @@
   #endif
 
   #include <memory>
-  #include <limits>
   #include <cstdlib>
   #include <cmath>
   #include <vector>
@@ -33,7 +33,7 @@ namespace quest
 namespace internal
 {
 template <typename T, int DIM>
-inline axom::StackArray<T, DIM> makeStackArray(T v = std::numeric_limits<T>::max())
+inline axom::StackArray<T, DIM> makeStackArray(T v = axom::numeric_limits<T>::max())
 {
   axom::StackArray<T, DIM> rval;
   for(int d = 0; d < DIM; ++d)
@@ -818,7 +818,7 @@ class MeshViewUtil
   MdIndices conduitIndicesToStackArray(
     const conduit::Node& node,
     const std::string& path,
-    axom::IndexType defaultVal = std::numeric_limits<axom::IndexType>::max()) const
+    axom::IndexType defaultVal = axom::numeric_limits<axom::IndexType>::max()) const
   {
     if(node.has_path(path))
     {
diff --git a/src/axom/quest/ScatteredInterpolation.hpp b/src/axom/quest/ScatteredInterpolation.hpp
index 27d2ecd974..4628bcd7b5 100644
--- a/src/axom/quest/ScatteredInterpolation.hpp
+++ b/src/axom/quest/ScatteredInterpolation.hpp
@@ -7,6 +7,7 @@
 #define QUEST_SCATTERED_INTERPOLATION_H_
 
 #include "axom/core.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/slic.hpp"
 #include "axom/sidre.hpp"
 #include "axom/spin.hpp"
@@ -19,7 +20,6 @@
 #include "conduit_blueprint.hpp"
 
 #include <cmath>
-#include <limits>
 
 namespace
 {
@@ -500,7 +500,7 @@ class ScatteredInterpolation
     conduit::Node& input_mesh,
     const std::string& input_field_name,
     const std::string& output_field_name,
-    const double INVALID_VALUE = std::numeric_limits<double>::quiet_NaN())
+    const double INVALID_VALUE = axom::numeric_limits<double>::quiet_NaN())
   {
     constexpr auto INVALID_INDEX = DelaunayTriangulation::INVALID_INDEX;
 
diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp
index 13d157db6c..dfe9793e99 100644
--- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp
+++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp
@@ -8,10 +8,11 @@
 
 #include "axom/config.hpp"
 #include "axom/core.hpp"
+#include "axom/core/NumericLimits.hpp"
+#include "axom/core/execution/runtime_policy.hpp"
 #include "axom/slic.hpp"
 #include "axom/primal.hpp"
 #include "axom/spin.hpp"
-#include "axom/core/execution/runtime_policy.hpp"
 
 #include "axom/fmt.hpp"
 
@@ -22,7 +23,6 @@
 #include "conduit_relay_io.hpp"
 
 #include <memory>
-#include <limits>
 #include <cstdlib>
 #include <cmath>
 #include <list>
@@ -254,7 +254,7 @@ class DistributedClosestPointImpl
     , m_mpiComm(MPI_COMM_NULL)
     , m_rank(-1)
     , m_nranks(-1)
-    , m_sqDistanceThreshold(std::numeric_limits<double>::max())
+    , m_sqDistanceThreshold(axom::numeric_limits<double>::max())
   { }
 
   virtual ~DistributedClosestPointImpl() { }
@@ -1042,9 +1042,9 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl
         cp_rank.fill(-1);
         cp_idx.fill(-1);
         cp_domidx.fill(-1);
-        const PointType nowhere(std::numeric_limits<double>::signaling_NaN());
+        const PointType nowhere(axom::numeric_limits<double>::signaling_NaN());
         cp_pos.fill(nowhere);
-        cp_dist.fill(std::numeric_limits<double>::signaling_NaN());
+        cp_dist.fill(axom::numeric_limits<double>::signaling_NaN());
       }
       auto query_inds = cp_idx.view();
       auto query_doms = cp_domidx.view();
diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp
index 935ce7b393..eed57dfddc 100644
--- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp
+++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp
@@ -11,6 +11,7 @@
 // Axom includes
 #include "axom/config.hpp"
 #include "axom/core.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/slic.hpp"
 #include "axom/primal.hpp"
 #include "axom/sidre.hpp"
@@ -35,7 +36,6 @@
 
 // C/C++ includes
 #include <string>
-#include <limits>
 #include <map>
 #include <vector>
 #include <cmath>
@@ -79,7 +79,7 @@ struct Input
 
   RuntimePolicy policy {RuntimePolicy::seq};
 
-  double distThreshold {std::numeric_limits<double>::max()};
+  double distThreshold {axom::numeric_limits<double>::max()};
 
   bool checkResults {false};
 
@@ -1195,8 +1195,8 @@ void computeDistancesAndDirections(BlueprintParticleMesh& queryMesh,
   using PointType = primal::Point<double, DIM>;
   using IndexSet = slam::PositionSet<>;
 
-  PointType nowhere(std::numeric_limits<double>::signaling_NaN());
-  const double nodist = std::numeric_limits<double>::signaling_NaN();
+  PointType nowhere(axom::numeric_limits<double>::signaling_NaN());
+  const double nodist = axom::numeric_limits<double>::signaling_NaN();
 
   queryMesh.registerNodalScalarField<double>(distanceField);
   queryMesh.registerNodalVectorField<double>(directionField);
diff --git a/src/axom/quest/examples/quest_marching_cubes_example.cpp b/src/axom/quest/examples/quest_marching_cubes_example.cpp
index 1cfac4f134..3527c5e90a 100644
--- a/src/axom/quest/examples/quest_marching_cubes_example.cpp
+++ b/src/axom/quest/examples/quest_marching_cubes_example.cpp
@@ -23,6 +23,7 @@
 
 // Axom includes
 #include "axom/core.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/slic.hpp"
 #include "axom/primal.hpp"
 #include "axom/mint/mesh/UnstructuredMesh.hpp"
@@ -49,7 +50,6 @@
 
 // C/C++ includes
 #include <string>
-#include <limits>
 #include <map>
 #include <vector>
 #include <cmath>
@@ -1479,8 +1479,8 @@ struct ContourTestBase
             cellMDMapper.toMultiIndex(parentCellId);
 
           // Compute min and max function values in the cell.
-          double minFcnValue = std::numeric_limits<double>::max();
-          double maxFcnValue = std::numeric_limits<double>::min();
+          double minFcnValue = axom::numeric_limits<double>::max();
+          double maxFcnValue = axom::numeric_limits<double>::min();
           constexpr short int cornerCount =
             (1 << DIM);  // Number of nodes in a cell.
           for(short int cornerId = 0; cornerId < cornerCount; ++cornerId)
diff --git a/src/axom/quest/interface/internal/QuestHelpers.cpp b/src/axom/quest/interface/internal/QuestHelpers.cpp
index 33b8317bd3..0ba10465c3 100644
--- a/src/axom/quest/interface/internal/QuestHelpers.cpp
+++ b/src/axom/quest/interface/internal/QuestHelpers.cpp
@@ -6,6 +6,7 @@
 #include "axom/quest/interface/internal/QuestHelpers.hpp"
 
 #include "axom/core.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/mint/mesh/UnstructuredMesh.hpp"
 
 // Quest includes
@@ -519,8 +520,8 @@ void compute_mesh_bounds(const mint::Mesh* mesh, double* lo, double* hi)
   // STEP 0: initialize lo,hi
   for(int i = 0; i < ndims; ++i)
   {
-    lo[i] = std::numeric_limits<double>::max();
-    hi[i] = std::numeric_limits<double>::lowest();
+    lo[i] = axom::numeric_limits<double>::max();
+    hi[i] = axom::numeric_limits<double>::lowest();
   }  // END for all dimensions
 
   // STEP 1: compute lo,hi
diff --git a/src/axom/quest/tests/quest_point_in_cell_mfem.cpp b/src/axom/quest/tests/quest_point_in_cell_mfem.cpp
index fff5f72c5a..89128df774 100644
--- a/src/axom/quest/tests/quest_point_in_cell_mfem.cpp
+++ b/src/axom/quest/tests/quest_point_in_cell_mfem.cpp
@@ -14,6 +14,7 @@
 
 #include "axom/config.hpp"
 #include "axom/core.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/mint.hpp"
 #include "axom/primal.hpp"
 #include "axom/spin.hpp"
@@ -195,7 +196,7 @@ class PointInCellTest : public ::testing::Test
 
     // compute minimal local mesh size
     mfem::Vector h0(fespace->GetNDofs());
-    h0 = std::numeric_limits<double>::infinity();
+    h0 = axom::numeric_limits<double>::infinity();
     {
       mfem::Array<int> dofs;
       for(int i = 0; i < fespace->GetNE(); i++)
diff --git a/src/axom/quest/tests/quest_pro_e_reader.cpp b/src/axom/quest/tests/quest_pro_e_reader.cpp
index 52f21f00d9..3c1de1804c 100644
--- a/src/axom/quest/tests/quest_pro_e_reader.cpp
+++ b/src/axom/quest/tests/quest_pro_e_reader.cpp
@@ -4,6 +4,7 @@
 // SPDX-License-Identifier: (BSD-3-Clause)
 
 #include "axom/core/utilities/FileUtilities.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/mint/utils/vtk_utils.hpp"  // for write_vtk
 #include "axom/quest/readers/ProEReader.hpp"
 #include "axom/slic.hpp"
@@ -132,13 +133,13 @@ TEST(quest_pro_e_reader, read_pro_e)
   {
     EXPECT_NEAR(x[inode],
                 x_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(y[inode],
                 y_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(z[inode],
                 z_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
   }  // END for all nodes
 
   // STEP 4: remove temporary Pro/E file
@@ -186,13 +187,13 @@ TEST(quest_pro_e_reader, read_pro_e_invbbox)
   {
     EXPECT_NEAR(x[inode],
                 x_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(y[inode],
                 y_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(z[inode],
                 z_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
   }  // END for all nodes
 
   // STEP 4: remove temporary Pro/E file
@@ -242,13 +243,13 @@ TEST(quest_pro_e_reader, read_pro_e_bbox_all)
   {
     EXPECT_NEAR(x[inode],
                 x_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(y[inode],
                 y_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(z[inode],
                 z_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
   }  // END for all nodes
 
   // STEP 4: remove temporary Pro/E file
@@ -299,13 +300,13 @@ TEST(quest_pro_e_reader, read_pro_e_bbox_some)
   {
     EXPECT_NEAR(x[inode],
                 x_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(y[inode],
                 y_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(z[inode],
                 z_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
   }  // END for all nodes
 
   // STEP 4: remove temporary Pro/E file
@@ -356,13 +357,13 @@ TEST(quest_pro_e_reader, read_pro_e_bbox_some_incl)
   {
     EXPECT_NEAR(x[inode],
                 x_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(y[inode],
                 y_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(z[inode],
                 z_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
   }  // END for all nodes
 
   // STEP 4: remove temporary Pro/E file
@@ -415,13 +416,13 @@ TEST(quest_pro_e_reader, read_pro_e_external)
   {
     EXPECT_NEAR(x[inode],
                 x_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(y[inode],
                 y_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(z[inode],
                 z_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
   }  // END for all nodes
 
   // STEP 4: remove temporary Pro/E file
@@ -434,7 +435,7 @@ TEST(quest_pro_e_reader, cup_pro_e)
 {
   constexpr int NUM_NODES = 171;
   constexpr int NUM_TETS = 574;
-  constexpr double EPS = std::numeric_limits<double>::epsilon();
+  constexpr double EPS = axom::numeric_limits<double>::epsilon();
 
   // STEP 0: Get Pro/E cup example file for testing
   namespace fs = axom::utilities::filesystem;
diff --git a/src/axom/quest/tests/quest_pro_e_reader_parallel.cpp b/src/axom/quest/tests/quest_pro_e_reader_parallel.cpp
index fcc9f1912b..e6a8c441fe 100644
--- a/src/axom/quest/tests/quest_pro_e_reader_parallel.cpp
+++ b/src/axom/quest/tests/quest_pro_e_reader_parallel.cpp
@@ -5,6 +5,8 @@
 
 #include "axom/config.hpp"
 
+#include "axom/core/NumericLimits.hpp"
+
 #include "axom/slic.hpp"
 
 #include "axom/quest/readers/PProEReader.hpp"
@@ -106,13 +108,13 @@ TEST(quest_pro_e_reader_parallel, read_file)
   {
     EXPECT_NEAR(x[inode],
                 x_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(y[inode],
                 y_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(z[inode],
                 z_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
   }  // END for all nodes
 
   // // STEP 4: remove temporary Pro/E file
@@ -173,13 +175,13 @@ TEST(quest_pro_e_reader_parallel, read_file_bbox)
   {
     EXPECT_NEAR(x[inode],
                 x_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(y[inode],
                 y_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(z[inode],
                 z_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
   }  // END for all nodes
 
   // // STEP 4: remove temporary Pro/E file
@@ -240,13 +242,13 @@ TEST(quest_pro_e_reader_parallel, read_file_bbox_incl)
   {
     EXPECT_NEAR(x[inode],
                 x_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(y[inode],
                 y_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(z[inode],
                 z_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
   }  // END for all nodes
 
   // // STEP 4: remove temporary Pro/E file
diff --git a/src/axom/quest/tests/quest_signed_distance.cpp b/src/axom/quest/tests/quest_signed_distance.cpp
index 2e10742e7e..484cfa1f84 100644
--- a/src/axom/quest/tests/quest_signed_distance.cpp
+++ b/src/axom/quest/tests/quest_signed_distance.cpp
@@ -4,6 +4,7 @@
 // SPDX-License-Identifier: (BSD-3-Clause)
 
 #include "axom/config.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/slic.hpp"
 #include "axom/mint.hpp"
 #include "axom/primal.hpp"
@@ -128,7 +129,7 @@ TEST(quest_signed_distance, sphere_test)
 
   double l1norm = 0.0;
   double l2norm = 0.0;
-  double linf = std::numeric_limits<double>::min();
+  double linf = axom::numeric_limits<double>::min();
 
   for(int inode = 0; inode < nnodes; ++inode)
   {
@@ -236,7 +237,7 @@ TEST(quest_signed_distance, sphere_test_with_normals)
 
   double l1norm = 0.0;
   double l2norm = 0.0;
-  double linf = std::numeric_limits<double>::min();
+  double linf = axom::numeric_limits<double>::min();
 
   for(int inode = 0; inode < nnodes; ++inode)
   {
@@ -369,7 +370,7 @@ void run_vectorized_sphere_test()
 
   double l1norm = 0.0;
   double l2norm = 0.0;
-  double linf = std::numeric_limits<double>::min();
+  double linf = axom::numeric_limits<double>::min();
 
   axom::Array<PointType> queryPts =
     axom::Array<PointType>(nnodes, nnodes, host_allocator);
@@ -552,7 +553,7 @@ TEST(quest_signed_distance, sphere_vec_device_custom_alloc)
 
   double l1norm = 0.0;
   double l2norm = 0.0;
-  double linf = std::numeric_limits<double>::min();
+  double linf = axom::numeric_limits<double>::min();
 
   axom::Array<PointType> queryPts =
     axom::Array<PointType>(nnodes, nnodes, host_allocator);
diff --git a/src/axom/quest/tests/quest_signed_distance_interface.cpp b/src/axom/quest/tests/quest_signed_distance_interface.cpp
index b94851a0c9..0efe5b2598 100644
--- a/src/axom/quest/tests/quest_signed_distance_interface.cpp
+++ b/src/axom/quest/tests/quest_signed_distance_interface.cpp
@@ -6,6 +6,7 @@
 // Axom utils
 #include "axom/config.hpp"
 #include "axom/core.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/mint.hpp"
 #include "axom/primal.hpp"
 #include "axom/slic.hpp"
@@ -480,7 +481,7 @@ TEST(quest_signed_distance_interface, analytic_sphere)
   // STEP 4: Compute signed distance
   double l1norm = 0.0;
   double l2norm = 0.0;
-  double linf = std::numeric_limits<double>::min();
+  double linf = axom::numeric_limits<double>::min();
   axom::IndexType nnodes = umesh->getNumberOfNodes();
   for(axom::IndexType inode = 0; inode < nnodes; ++inode)
   {
@@ -577,7 +578,7 @@ TEST(quest_signed_distance_interface, analytic_sphere_with_closest_pt_and_normal
   // STEP 4: Compute signed distance
   double l1norm = 0.0;
   double l2norm = 0.0;
-  double linf = std::numeric_limits<double>::min();
+  double linf = axom::numeric_limits<double>::min();
   axom::IndexType nnodes = umesh->getNumberOfNodes();
   for(axom::IndexType inode = 0; inode < nnodes; ++inode)
   {
diff --git a/src/axom/quest/tests/quest_stl_reader.cpp b/src/axom/quest/tests/quest_stl_reader.cpp
index 071ad2805c..b483463d44 100644
--- a/src/axom/quest/tests/quest_stl_reader.cpp
+++ b/src/axom/quest/tests/quest_stl_reader.cpp
@@ -6,6 +6,7 @@
 #include "axom/quest/readers/STLReader.hpp"
 #include "axom/mint/mesh/UnstructuredMesh.hpp"
 #include "axom/slic.hpp"
+#include "axom/core/NumericLimits.hpp"
 
 // gtest includes
 #include "gtest/gtest.h"
@@ -14,7 +15,6 @@
 #include <cstdio>
 #include <string>
 #include <fstream>
-#include <limits>
 
 // namespace aliases
 namespace mint = axom::mint;
@@ -133,13 +133,13 @@ TEST(quest_stl_reader, read_stl)
   {
     EXPECT_NEAR(x[inode],
                 x_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(y[inode],
                 y_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(z[inode],
                 z_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
   }  // END for all nodes
 
   // STEP 4: remove temporary STL file
@@ -192,13 +192,13 @@ TEST(quest_stl_reader, read_stl_external)
   {
     EXPECT_NEAR(x[inode],
                 x_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(y[inode],
                 y_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
     EXPECT_NEAR(z[inode],
                 z_expected[inode],
-                std::numeric_limits<double>::epsilon());
+                axom::numeric_limits<double>::epsilon());
   }  // END for all nodes
 
   // STEP 4: remove temporary STL file
diff --git a/src/axom/sidre/tests/sidre_types_C.cpp b/src/axom/sidre/tests/sidre_types_C.cpp
index 2c8c90ce2b..6ed8f43066 100644
--- a/src/axom/sidre/tests/sidre_types_C.cpp
+++ b/src/axom/sidre/tests/sidre_types_C.cpp
@@ -5,9 +5,8 @@
 
 #include "gtest/gtest.h"
 
-#include <limits>
-
 #include "axom/core/Types.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/sidre.hpp"
 
 using axom::sidre::DataType;
@@ -24,21 +23,21 @@ template <typename CommonType, typename SidreType>
 void testTypesForEquality()
 {
   // check if both are integral types
-  bool com_is_int = std::numeric_limits<CommonType>::is_integer;
-  bool sid_is_int = std::numeric_limits<SidreType>::is_integer;
+  bool com_is_int = axom::numeric_limits<CommonType>::is_integer;
+  bool sid_is_int = axom::numeric_limits<SidreType>::is_integer;
   EXPECT_EQ(com_is_int, sid_is_int);
 
   // check if both are signed (or both are unsigned)
-  bool com_is_signed = std::numeric_limits<CommonType>::is_signed;
-  bool sid_is_signed = std::numeric_limits<SidreType>::is_signed;
+  bool com_is_signed = axom::numeric_limits<CommonType>::is_signed;
+  bool sid_is_signed = axom::numeric_limits<SidreType>::is_signed;
   EXPECT_EQ(com_is_signed, sid_is_signed);
 
   // check that both have same number of bytes
   EXPECT_EQ(sizeof(CommonType), sizeof(SidreType));
 
   // check that both have same number of digits
-  EXPECT_EQ(std::numeric_limits<CommonType>::digits,
-            std::numeric_limits<SidreType>::digits);
+  EXPECT_EQ(axom::numeric_limits<CommonType>::digits,
+            axom::numeric_limits<SidreType>::digits);
 }
 
 }  // namespace
diff --git a/src/axom/slam/policies/IndirectionPolicies.hpp b/src/axom/slam/policies/IndirectionPolicies.hpp
index 8b2bceda94..66675de942 100644
--- a/src/axom/slam/policies/IndirectionPolicies.hpp
+++ b/src/axom/slam/policies/IndirectionPolicies.hpp
@@ -34,6 +34,7 @@
 
 #include "axom/core/Macros.hpp"
 #include "axom/core/Array.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/slic/interface/slic.hpp"
 
 namespace axom
@@ -303,7 +304,7 @@ struct CArrayIndirectionBase
 
   constexpr PositionType size() const
   {
-    return std::numeric_limits<PositionType>::max();
+    return axom::numeric_limits<PositionType>::max();
   }
 
 private:
diff --git a/src/axom/spin/MortonIndex.hpp b/src/axom/spin/MortonIndex.hpp
index ce9d8ad5ef..e23c98a71d 100644
--- a/src/axom/spin/MortonIndex.hpp
+++ b/src/axom/spin/MortonIndex.hpp
@@ -20,10 +20,10 @@
 #include "axom/config.hpp"
 #include "axom/core/Types.hpp"
 #include "axom/core/Macros.hpp"  // defines AXOM_STATIC_ASSERT
+#include "axom/core/NumericLimits.hpp"
 #include "axom/primal/geometry/Point.hpp"
 
 #include <type_traits>
-#include <limits>  // for numeric_limits
 
 namespace
 {
@@ -257,10 +257,10 @@ struct Mortonizer<CoordType, MortonIndexType, 2>
     NDIM = 2,
 
     /*! The number of bits in a CoordType  */
-    COORD_BITS = std::numeric_limits<CoordType>::digits,
+    COORD_BITS = axom::numeric_limits<CoordType>::digits,
 
     /*! The number of bits in a MortonIndex  */
-    MORTON_BITS = std::numeric_limits<MortonIndexType>::digits,
+    MORTON_BITS = axom::numeric_limits<MortonIndexType>::digits,
 
     /*! The number of representable Morton bits per dimension */
     MB_PER_DIM = MORTON_BITS / NDIM,
@@ -399,10 +399,10 @@ struct Mortonizer<CoordType, MortonIndexType, 3>
     NDIM = 3,
 
     /*! The number of bits in a CoordType  */
-    COORD_BITS = std::numeric_limits<CoordType>::digits,
+    COORD_BITS = axom::numeric_limits<CoordType>::digits,
 
     /*! The number of bits in a MortonIndex  */
-    MORTON_BITS = std::numeric_limits<MortonIndexType>::digits,
+    MORTON_BITS = axom::numeric_limits<MortonIndexType>::digits,
 
     /*! The number of representable morton bits per dimension */
     MB_PER_DIM = MORTON_BITS / NDIM,
diff --git a/src/axom/spin/OctreeBase.hpp b/src/axom/spin/OctreeBase.hpp
index 3bf0e6ac16..9aae342063 100644
--- a/src/axom/spin/OctreeBase.hpp
+++ b/src/axom/spin/OctreeBase.hpp
@@ -12,6 +12,7 @@
 #define AXOM_SPIN_OCTREE_BASE__HPP_
 
 #include "axom/config.hpp"
+#include "axom/core/NumericLimits.hpp"
 #include "axom/slic.hpp"
 #include "axom/slam.hpp"
 #include "axom/primal.hpp"
@@ -140,7 +141,7 @@ class OctreeBase
   using GridVec = primal::Vector<CoordType, DIM>;
 
   using MAX_LEVEL_SIZE =
-    slam::policies::CompileTimeSize<CoordType, std::numeric_limits<CoordType>::digits>;
+    slam::policies::CompileTimeSize<CoordType, axom::numeric_limits<CoordType>::digits>;
   using OctreeLevels = slam::OrderedSet<CoordType, CoordType, MAX_LEVEL_SIZE>;
 
   using OctreeLevelType = OctreeLevel<DIM, BlockDataType>;
diff --git a/src/axom/spin/SparseOctreeLevel.hpp b/src/axom/spin/SparseOctreeLevel.hpp
index 4abe500042..5756f2b3ee 100644
--- a/src/axom/spin/SparseOctreeLevel.hpp
+++ b/src/axom/spin/SparseOctreeLevel.hpp
@@ -9,6 +9,7 @@
 #include "axom/config.hpp"
 
 #include "axom/core.hpp"
+#include "axom/core/NumericLimits.hpp"
 
 #include "axom/primal/geometry/Point.hpp"
 
@@ -72,7 +73,7 @@ struct BroodRepresentationTraits
   {
 #if defined(AXOM_USE_SPARSEHASH)
     const PointRepresenationType maxVal =
-      std::numeric_limits<PointRepresenationType>::max();
+      axom::numeric_limits<PointRepresenationType>::max();
     map.set_empty_key(maxVal);
     map.set_deleted_key(maxVal - 1);
 #else
@@ -125,7 +126,7 @@ struct BroodRepresentationTraits<CoordType, DIM, BroodDataType, primal::Point<Co
   static void initializeMap(MapType& map)
   {
 #if defined(AXOM_USE_SPARSEHASH)
-    CoordType maxCoord = std::numeric_limits<CoordType>::max();
+    CoordType maxCoord = axom::numeric_limits<CoordType>::max();
     GridPt maxPt(maxCoord);
     map.set_empty_key(maxPt);
 
diff --git a/src/axom/spin/UniformGrid.hpp b/src/axom/spin/UniformGrid.hpp
index 30c495ba58..39ec3211f9 100644
--- a/src/axom/spin/UniformGrid.hpp
+++ b/src/axom/spin/UniformGrid.hpp
@@ -9,6 +9,7 @@
 #include "axom/core/utilities/Utilities.hpp"
 #include "axom/core/execution/for_all.hpp"
 #include "axom/core/Array.hpp"
+#include "axom/core/NumericLimits.hpp"
 
 #include "axom/slic/interface/slic.hpp"
 
@@ -475,8 +476,8 @@ void UniformGrid<T, NDIMS, ExecSpace, StoragePolicy>::initialize(
   // get the global bounding box of all the objects
 #ifdef AXOM_USE_RAJA
   using reduce_pol = typename axom::execution_space<ExecSpace>::reduce_policy;
-  double infinity = std::numeric_limits<double>::max();
-  double neg_infinity = std::numeric_limits<double>::lowest();
+  double infinity = axom::numeric_limits<double>::max();
+  double neg_infinity = axom::numeric_limits<double>::lowest();
 
   using ReduceMin = RAJA::ReduceMin<reduce_pol, double>;
   using ReduceMax = RAJA::ReduceMax<reduce_pol, double>;
diff --git a/src/axom/spin/internal/linear_bvh/build_radix_tree.hpp b/src/axom/spin/internal/linear_bvh/build_radix_tree.hpp
index 9937b8ed73..b86d606fad 100644
--- a/src/axom/spin/internal/linear_bvh/build_radix_tree.hpp
+++ b/src/axom/spin/internal/linear_bvh/build_radix_tree.hpp
@@ -13,6 +13,7 @@
 #include "axom/core/AnnotationMacros.hpp"
 #include "axom/core/utilities/Utilities.hpp"
 #include "axom/core/utilities/BitUtilities.hpp"
+#include "axom/core/NumericLimits.hpp"
 
 #include "axom/slic/interface/slic.hpp"
 
@@ -117,8 +118,8 @@ primal::BoundingBox<FloatType, NDIMS> reduce(
 
   primal::Point<FloatType, NDIMS> min_pt, max_pt;
 
-  FloatType infinity = std::numeric_limits<FloatType>::max();
-  FloatType neg_infinity = std::numeric_limits<FloatType>::lowest();
+  FloatType infinity = axom::numeric_limits<FloatType>::max();
+  FloatType neg_infinity = axom::numeric_limits<FloatType>::lowest();
 
   for(int dim = 0; dim < NDIMS; dim++)
   {
diff --git a/src/axom/spin/tests/spin_morton.cpp b/src/axom/spin/tests/spin_morton.cpp
index 1f1ed34d98..a82a0edeab 100644
--- a/src/axom/spin/tests/spin_morton.cpp
+++ b/src/axom/spin/tests/spin_morton.cpp
@@ -5,13 +5,14 @@
 
 #include "gtest/gtest.h"
 
+#include "axom/core/NumericLimits.hpp"
+
 #include "axom/spin/MortonIndex.hpp"
 
 #include "axom/primal/geometry/Point.hpp"
 #include "axom/slic.hpp"
 
 #include <cstdlib>
-#include <limits>
 
 // Uncomment the line below for true randomized points
 #ifndef MORTON_TESTER_SHOULD_SEED
@@ -36,7 +37,7 @@ CoordType randomInt(CoordType beg, CoordType end)
 
   if(range == 0)
   {
-    range = std::numeric_limits<CoordType>::max();
+    range = axom::numeric_limits<CoordType>::max();
   }
 
   return (std::rand() % range) + beg;
@@ -65,7 +66,7 @@ TEST(spin_morton, test_max_set_bit)
   axom::spin::Mortonizer<CoordType, MortonIndexType, 2> morton2;
   EXPECT_EQ(morton2.maxSetBit(0), 0);
 
-  int maxBit = std::numeric_limits<CoordType>::digits;
+  int maxBit = axom::numeric_limits<CoordType>::digits;
   for(int i = 0; i <= maxBit; ++i)
   {
     int val = 1 << i;
@@ -135,7 +136,7 @@ void testMortonizer()
   int maxBits =
     axom::spin::Mortonizer<CoordType, MortonIndexType, DIM>::maxBitsPerCoord();
   SLIC_INFO(
-    "\tMax bits per dimension: " << std::numeric_limits<CoordType>::digits);
+    "\tMax bits per dimension: " << axom::numeric_limits<CoordType>::digits);
   SLIC_INFO("\tMax unique bits per dimension: " << maxBits);
 
   int maxIter = std::min(1 << (maxBits - 1), MAX_ITER);