diff --git a/.docs/Notebooks/sfrpackage_example.py b/.docs/Notebooks/sfrpackage_example.py
index 292f7d6fec..7aa328a13d 100644
--- a/.docs/Notebooks/sfrpackage_example.py
+++ b/.docs/Notebooks/sfrpackage_example.py
@@ -204,7 +204,6 @@
     raise ValueError("Failed to run.")
 
 # ### Load SFR formated water balance output into pandas dataframe using the `SfrFile` class
-# * requires the **pandas** library
 
 sfr_outfile = os.path.join(
     "..", "..", "examples", "data", "sfr_examples", "test1ss.flw"
diff --git a/autotest/test_export.py b/autotest/test_export.py
index a21c060c11..9f4eac4d32 100644
--- a/autotest/test_export.py
+++ b/autotest/test_export.py
@@ -160,7 +160,7 @@ def test_output_helper_shapefile_export(
     )
 
 
-@requires_pkg("pandas", "shapefile")
+@requires_pkg("shapefile")
 @pytest.mark.slow
 def test_freyberg_export(function_tmpdir, example_data_path):
     # steady state
@@ -254,7 +254,7 @@ def test_freyberg_export(function_tmpdir, example_data_path):
                 assert part.read_text() == wkt
 
 
-@requires_pkg("pandas", "shapefile")
+@requires_pkg("shapefile")
 @pytest.mark.parametrize("missing_arrays", [True, False])
 @pytest.mark.slow
 def test_disu_export(function_tmpdir, missing_arrays):
@@ -485,7 +485,7 @@ def test_shapefile_ibound(function_tmpdir, example_data_path):
     shape.close()
 
 
-@requires_pkg("pandas", "shapefile")
+@requires_pkg("shapefile")
 @pytest.mark.slow
 @pytest.mark.parametrize("namfile", namfiles())
 def test_shapefile(function_tmpdir, namfile):
@@ -510,7 +510,7 @@ def test_shapefile(function_tmpdir, namfile):
     ), f"wrong number of records in shapefile {fnc_name}"
 
 
-@requires_pkg("pandas", "shapefile")
+@requires_pkg("shapefile")
 @pytest.mark.slow
 @pytest.mark.parametrize("namfile", namfiles())
 def test_shapefile_export_modelgrid_override(function_tmpdir, namfile):
@@ -1446,7 +1446,7 @@ def test_vtk_vertex(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf2005")
-@requires_pkg("pandas", "vtk")
+@requires_pkg("vtk")
 def test_vtk_pathline(function_tmpdir, example_data_path):
     from vtkmodules.vtkIOLegacy import vtkUnstructuredGridReader
 
diff --git a/autotest/test_grid.py b/autotest/test_grid.py
index 246f9e48b9..c3dafe5557 100644
--- a/autotest/test_grid.py
+++ b/autotest/test_grid.py
@@ -1,5 +1,5 @@
-import re
 import os
+import re
 import warnings
 from warnings import warn
 
@@ -22,7 +22,6 @@
 from flopy.utils.triangle import Triangle
 from flopy.utils.voronoi import VoronoiGrid
 
-
 HAS_PYPROJ = has_pkg("pyproj")
 if HAS_PYPROJ:
     import pyproj
diff --git a/autotest/test_hydmodfile.py b/autotest/test_hydmodfile.py
index 8e459a8241..76b70f1bb2 100644
--- a/autotest/test_hydmodfile.py
+++ b/autotest/test_hydmodfile.py
@@ -1,6 +1,7 @@
 import os
 
 import numpy as np
+import pandas as pd
 import pytest
 from modflow_devtools.markers import requires_pkg
 from modflow_devtools.misc import has_pkg
@@ -118,31 +119,22 @@ def test_hydmodfile_read(hydmod_model_path):
         len(data.dtype.names) == nitems + 1
     ), f"data column length is not {len(nitems + 1)}"
 
-    if has_pkg("pandas"):
-        import pandas as pd
-
-        for idx in range(ntimes):
-            df = h.get_dataframe(idx=idx, timeunit="S")
-            assert isinstance(df, pd.DataFrame), "A DataFrame was not returned"
-            assert df.shape == (1, 9), "data shape is not (1, 9)"
-
-        for time in times:
-            df = h.get_dataframe(totim=time, timeunit="S")
-            assert isinstance(df, pd.DataFrame), "A DataFrame was not returned"
-            assert df.shape == (1, 9), "data shape is not (1, 9)"
+    for idx in range(ntimes):
+        df = h.get_dataframe(idx=idx, timeunit="S")
+        assert isinstance(df, pd.DataFrame), "A DataFrame was not returned"
+        assert df.shape == (1, 9), "data shape is not (1, 9)"
 
-        df = h.get_dataframe(timeunit="S")
+    for time in times:
+        df = h.get_dataframe(totim=time, timeunit="S")
         assert isinstance(df, pd.DataFrame), "A DataFrame was not returned"
-        assert df.shape == (101, 9), "data shape is not (101, 9)"
-    else:
-        print("pandas not available...")
-        pass
+        assert df.shape == (1, 9), "data shape is not (1, 9)"
 
+    df = h.get_dataframe(timeunit="S")
+    assert isinstance(df, pd.DataFrame), "A DataFrame was not returned"
+    assert df.shape == (101, 9), "data shape is not (101, 9)"
 
-@requires_pkg("pandas")
-def test_mf6obsfile_read(mf6_obs_model_path):
-    import pandas as pd
 
+def test_mf6obsfile_read(mf6_obs_model_path):
     txt = "binary mf6 obs"
     files = ["maw_obs.gitbin", "maw_obs.gitcsv"]
     binfile = [True, False]
diff --git a/autotest/test_listbudget.py b/autotest/test_listbudget.py
index ffc4a094a3..7f39e8d089 100644
--- a/autotest/test_listbudget.py
+++ b/autotest/test_listbudget.py
@@ -2,6 +2,7 @@
 import warnings
 
 import numpy as np
+import pandas as pd
 import pytest
 from modflow_devtools.markers import requires_pkg
 from modflow_devtools.misc import has_pkg
@@ -54,14 +55,9 @@ def test_mflistfile(example_data_path):
     cum = mflist.get_cumulative(names="PERCENT_DISCREPANCY")
     assert isinstance(cum, np.ndarray)
 
-    if not has_pkg("pandas"):
-        return
-
-    import pandas
-
     df_flx, df_vol = mflist.get_dataframes(start_datetime=None)
-    assert isinstance(df_flx, pandas.DataFrame)
-    assert isinstance(df_vol, pandas.DataFrame)
+    assert isinstance(df_flx, pd.DataFrame)
+    assert isinstance(df_vol, pd.DataFrame)
 
     # test get runtime
     runtime = mflist.get_model_runtime(units="hours")
@@ -115,10 +111,7 @@ def test_mflist_reducedpumping_fail(example_data_path):
         mflist.get_reduced_pumping()
 
 
-@requires_pkg("pandas")
 def test_mtlist(example_data_path):
-    import pandas as pd
-
     mt_dir = example_data_path / "mt3d_test"
     mt = MtListBudget(mt_dir / "mcomp.list")
     df_gw, df_sw = mt.parse(forgive=False, diff=False, start_datetime=None)
diff --git a/autotest/test_mnw.py b/autotest/test_mnw.py
index d2b4c0ea04..266c55a14b 100644
--- a/autotest/test_mnw.py
+++ b/autotest/test_mnw.py
@@ -2,6 +2,7 @@
 import shutil
 
 import numpy as np
+import pandas as pd
 import pytest
 from modflow_devtools.markers import requires_pkg
 
@@ -293,13 +294,11 @@ def test_make_package(function_tmpdir):
     )
 
 
-@requires_pkg("pandas")
 def test_mnw2_create_file(function_tmpdir):
     """
     Test for issue #556, Mnw2 crashed if wells have
     multiple node lengths
     """
-    import pandas as pd
 
     mf = Modflow("test_mfmnw2", exe_name="mf2005")
     ws = function_tmpdir
diff --git a/autotest/test_mp5.py b/autotest/test_mp5.py
index 9101d0809d..7c4515444b 100644
--- a/autotest/test_mp5.py
+++ b/autotest/test_mp5.py
@@ -1,6 +1,7 @@
 import os
 
 import numpy as np
+import pandas as pd
 from autotest.test_mp6 import eval_timeseries
 from matplotlib import pyplot as plt
 from modflow_devtools.markers import requires_pkg
@@ -10,7 +11,6 @@
 from flopy.utils import EndpointFile, PathlineFile
 
 
-@requires_pkg("pandas")
 def test_mp5_load(function_tmpdir, example_data_path):
     # load the base freyberg model
     freyberg_ws = example_data_path / "freyberg"
@@ -68,7 +68,6 @@ def test_mp5_load(function_tmpdir, example_data_path):
     plt.close()
 
 
-@requires_pkg("pandas")
 def test_mp5_timeseries_load(example_data_path):
     pth = str(example_data_path / "mp5")
     files = [
diff --git a/autotest/test_mp6.py b/autotest/test_mp6.py
index 4637876e01..6f37de0369 100644
--- a/autotest/test_mp6.py
+++ b/autotest/test_mp6.py
@@ -3,6 +3,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 import pytest
 from autotest.conftest import get_example_data_path
 from autotest.test_mp6_cases import Mp6Cases1, Mp6Cases2
@@ -110,12 +111,8 @@ def test_mpsim(function_tmpdir, mp6_test_path):
     )
     mp.write_input()
 
-    use_pandas_combs = [False]  # test StartingLocationsFile._write_wo_pandas
-    if has_pkg("pandas"):
-        # test StartingLocationsFile._write_particle_data_with_pandas
-        use_pandas_combs.append(True)
-
-    for use_pandas in use_pandas_combs:
+    # test StartingLocationsFile._write_wo_pandas
+    for use_pandas in [True, False]:
         sim = Modpath6Sim(model=mp)
         # starting locations file
         stl = StartingLocationsFile(model=mp, use_pandas=use_pandas)
@@ -135,7 +132,7 @@ def test_mpsim(function_tmpdir, mp6_test_path):
         assert stllines[6].strip().split()[-1] == "p2"
 
 
-@requires_pkg("pandas", "shapefile", "shapely")
+@requires_pkg("shapefile", "shapely")
 def test_get_destination_data(function_tmpdir, mp6_test_path):
     copy_modpath_files(mp6_test_path, function_tmpdir, "EXAMPLE.")
     copy_modpath_files(mp6_test_path, function_tmpdir, "EXAMPLE-3.")
@@ -307,7 +304,6 @@ def test_get_destination_data(function_tmpdir, mp6_test_path):
     pthobj.write_shapefile(shpname=fpth, direction="ending", mg=mg4)
 
 
-@requires_pkg("pandas")
 def test_loadtxt(function_tmpdir, mp6_test_path):
     copy_modpath_files(mp6_test_path, function_tmpdir, "EXAMPLE-3.")
 
@@ -324,7 +320,6 @@ def test_loadtxt(function_tmpdir, mp6_test_path):
 
 
 @requires_exe("mf2005")
-@requires_pkg("pandas")
 def test_modpath(function_tmpdir, example_data_path):
     pth = example_data_path / "freyberg"
     mfnam = "freyberg.nam"
@@ -480,7 +475,6 @@ def test_modpath(function_tmpdir, example_data_path):
         plt.close()
 
 
-@requires_pkg("pandas")
 def test_mp6_timeseries_load(example_data_path):
     pth = example_data_path / "mp5"
     files = [
diff --git a/autotest/test_mp7.py b/autotest/test_mp7.py
index b51edf624e..540b6b9c25 100644
--- a/autotest/test_mp7.py
+++ b/autotest/test_mp7.py
@@ -254,7 +254,6 @@ def test_default_modpath(ex01b_mf6_model):
 
 
 @requires_exe("mf6", "mp7")
-@requires_pkg("pandas")
 def test_faceparticles_is1(ex01b_mf6_model):
     sim, function_tmpdir = ex01b_mf6_model
 
@@ -441,7 +440,6 @@ def test_facenode_is2a(ex01b_mf6_model):
 
 
 @requires_exe("mf6", "mp7")
-@requires_pkg("pandas")
 def test_cellparticles_is1(ex01b_mf6_model):
     sim, function_tmpdir = ex01b_mf6_model
     grid = sim.get_model(ex01b_mf6_model_name).modelgrid
@@ -817,7 +815,6 @@ def test_pathline_output(function_tmpdir):
     assert maxid0 == maxid1, msg
 
 
-@requires_pkg("pandas")
 @requires_exe("mf2005", "mf6", "mp7")
 def test_endpoint_output(function_tmpdir):
     case_mf2005 = Mp7Cases.mp7_mf2005(function_tmpdir)
diff --git a/autotest/test_plot.py b/autotest/test_plot.py
index 085eabd693..04f320721a 100644
--- a/autotest/test_plot.py
+++ b/autotest/test_plot.py
@@ -386,7 +386,6 @@ def modpath_model(function_tmpdir, example_data_path):
     return ml, mp, sim
 
 
-@requires_pkg("pandas")
 @requires_exe("mf2005", "mp6")
 def test_xc_plot_particle_pathlines(modpath_model):
     ml, mp, sim = modpath_model
@@ -407,7 +406,6 @@ def test_xc_plot_particle_pathlines(modpath_model):
     assert len(pth._paths) == 6
 
 
-@requires_pkg("pandas")
 @requires_exe("mf2005", "mp6")
 def test_map_plot_particle_endpoints(modpath_model):
     ml, mp, sim = modpath_model
diff --git a/autotest/test_sfr.py b/autotest/test_sfr.py
index 19b09d0fb2..42630e1880 100644
--- a/autotest/test_sfr.py
+++ b/autotest/test_sfr.py
@@ -373,7 +373,7 @@ def test_const(sfr_data):
     assert True
 
 
-@requires_pkg("pandas", "shapefile", "shapely")
+@requires_pkg("shapefile", "shapely")
 def test_export(function_tmpdir, sfr_data):
     m = Modflow()
     dis = ModflowDis(m, 1, 10, 10, lenuni=2, itmuni=4)
@@ -666,7 +666,6 @@ def test_assign_layers(function_tmpdir):
 
 
 @requires_exe("mf2005")
-@requires_pkg("pandas")
 def test_SfrFile(function_tmpdir, sfr_examples_path, mf2005_model_path):
     common_names = [
         "layer",
@@ -693,13 +692,12 @@ def test_SfrFile(function_tmpdir, sfr_examples_path, mf2005_model_path):
         "gw_head",
     ], sfrout.names
     assert sfrout.times == [(0, 0), (49, 1)], sfrout.times
-    # will be None if pandas is not installed
-    if sfrout.pd is not None:
-        df = sfrout.get_dataframe()
-        assert df.layer.values[0] == 1
-        assert df.column.values[0] == 169
-        assert df.Cond.values[0] == 74510.0
-        assert df.gw_head.values[3] == 1.288e03
+
+    df = sfrout.get_dataframe()
+    assert df.layer.values[0] == 1
+    assert df.column.values[0] == 169
+    assert df.Cond.values[0] == 74510.0
+    assert df.gw_head.values[3] == 1.288e03
 
     sfrout = SfrFile(sfr_examples_path / "test1tr.flw")
     assert sfrout.ncol == 16, sfrout.ncol
@@ -737,10 +735,9 @@ def test_SfrFile(function_tmpdir, sfr_examples_path, mf2005_model_path):
         (49, 1),
     ]
     assert sfrout.times == expected_times, sfrout.times
-    if sfrout.pd is not None:
-        df = sfrout.get_dataframe()
-        assert df.gradient.values[-1] == 5.502e-02
-        assert df.shape == (1080, 20)
+    df = sfrout.get_dataframe()
+    assert df.gradient.values[-1] == 5.502e-02
+    assert df.shape == (1080, 20)
 
     ml = Modflow.load(
         "test1tr.nam", model_ws=mf2005_model_path, exe_name="mf2005"
diff --git a/autotest/test_str.py b/autotest/test_str.py
index b54149137a..85c23333f2 100644
--- a/autotest/test_str.py
+++ b/autotest/test_str.py
@@ -14,7 +14,6 @@
 
 
 @requires_exe("mf2005")
-@requires_pkg("pandas")
 def test_str_issue1164(function_tmpdir, example_data_path):
     mf2005_model_path = example_data_path / "mf2005_test"
     m = Modflow.load(
diff --git a/autotest/test_swr_binaryread.py b/autotest/test_swr_binaryread.py
index 228127a12b..02bf20a96b 100644
--- a/autotest/test_swr_binaryread.py
+++ b/autotest/test_swr_binaryread.py
@@ -1,4 +1,5 @@
 # Test SWR binary read functionality
+import pandas as pd
 import pytest
 from modflow_devtools.misc import has_pkg
 
@@ -446,21 +447,16 @@ def test_swr_binary_obs(swr_test_path, ipos):
         ), "SwrObs data does not have nobs + 1"
 
     # test get_dataframes()
-    if has_pkg("pandas"):
-        import pandas as pd
-
-        for idx in range(ntimes):
-            df = sobj.get_dataframe(idx=idx, timeunit="S")
-            assert isinstance(df, pd.DataFrame), "A DataFrame was not returned"
-            assert df.shape == (1, nobs + 1), "data shape is not (1, 10)"
-
-        for time in times:
-            df = sobj.get_dataframe(totim=time, timeunit="S")
-            assert isinstance(df, pd.DataFrame), "A DataFrame was not returned"
-            assert df.shape == (1, nobs + 1), "data shape is not (1, 10)"
+    for idx in range(ntimes):
+        df = sobj.get_dataframe(idx=idx, timeunit="S")
+        assert isinstance(df, pd.DataFrame), "A DataFrame was not returned"
+        assert df.shape == (1, nobs + 1), "data shape is not (1, 10)"
 
-        df = sobj.get_dataframe(timeunit="S")
+    for time in times:
+        df = sobj.get_dataframe(totim=time, timeunit="S")
         assert isinstance(df, pd.DataFrame), "A DataFrame was not returned"
-        assert df.shape == (336, nobs + 1), "data shape is not (336, 10)"
-    else:
-        print("pandas not available...")
+        assert df.shape == (1, nobs + 1), "data shape is not (1, 10)"
+
+    df = sobj.get_dataframe(timeunit="S")
+    assert isinstance(df, pd.DataFrame), "A DataFrame was not returned"
+    assert df.shape == (336, nobs + 1), "data shape is not (336, 10)"
diff --git a/autotest/test_util_2d_and_3d.py b/autotest/test_util_2d_and_3d.py
index ea262202fb..2080458033 100644
--- a/autotest/test_util_2d_and_3d.py
+++ b/autotest/test_util_2d_and_3d.py
@@ -1,6 +1,7 @@
 import os
 
 import numpy as np
+import pandas as pd
 import pytest
 from modflow_devtools.markers import requires_pkg
 
@@ -439,7 +440,6 @@ def test_append_mflist(function_tmpdir):
     ml.write_input()
 
 
-@requires_pkg("pandas")
 def test_mflist(function_tmpdir, example_data_path):
     model = Modflow(model_ws=function_tmpdir)
     dis = ModflowDis(model, 10, 10, 10, 10)
@@ -609,13 +609,10 @@ def test_util3d_reset():
     ml.bas6.strt = arr
 
 
-@requires_pkg("pandas")
 def test_mflist_fromfile(function_tmpdir):
     """test that when a file is passed to stress period data,
     the .array attribute will load the file
     """
-    import pandas as pd
-
     wel_data = pd.DataFrame(
         [(0, 1, 2, -50.0), (0, 5, 5, -50.0)], columns=["k", "i", "j", "flux"]
     )
diff --git a/autotest/test_zonbud_utility.py b/autotest/test_zonbud_utility.py
index 57b3a7d49f..4dd16c70c6 100644
--- a/autotest/test_zonbud_utility.py
+++ b/autotest/test_zonbud_utility.py
@@ -1,6 +1,7 @@
 import os
 
 import numpy as np
+import pandas as pd
 import pytest
 from modflow_devtools.markers import requires_exe, requires_pkg
 
@@ -207,7 +208,6 @@ def test_zonbud_readwrite_zbarray(function_tmpdir):
     assert np.array_equal(x, z), "Input and output arrays do not match."
 
 
-@requires_pkg("pandas")
 def test_dataframes(cbc_f, zon_f):
     zon = ZoneBudget.read_zone_file(zon_f)
     cmd = ZoneBudget(cbc_f, zon, totim=1095.0)
@@ -233,11 +233,8 @@ def test_get_model_shape(cbc_f, zon_f):
     ).get_model_shape()
 
 
-@requires_pkg("pandas")
 @pytest.mark.parametrize("rtol", [1e-2])
 def test_zonbud_active_areas_zone_zero(loadpth, cbc_f, rtol):
-    import pandas as pd
-
     # Read ZoneBudget executable output and reformat
     zbud_f = loadpth / "zonef_mlt_active_zone_0.2.csv"
     zbud = pd.read_csv(zbud_f)
@@ -284,10 +281,7 @@ def test_read_zone_file(function_tmpdir):
 
 @pytest.mark.mf6
 @requires_exe("mf6")
-@requires_pkg("pandas")
 def test_zonebudget_6(function_tmpdir, example_data_path):
-    import pandas as pd
-
     exe_name = "mf6"
     zb_exe_name = "zbud6"
 
diff --git a/docs/flopy_method_dependencies.md b/docs/flopy_method_dependencies.md
index 28ddaa89f3..72841add05 100644
--- a/docs/flopy_method_dependencies.md
+++ b/docs/flopy_method_dependencies.md
@@ -11,17 +11,6 @@ Additional dependencies to use optional FloPy helper methods are listed below.
 | `.interpolate()` in `flopy.utils.reference` `SpatialReference` class                 | **scipy.interpolate**                                                    |
 | `.interpolate()` in `flopy.mf6.utils.reference` `StructuredSpatialReference` class   | **scipy.interpolate**                                                    |
 | `._parse_units_from_proj4()` in `flopy.utils.reference` `SpatialReference` class     | **pyproj**                                                               |
-| `.get_dataframes()` in `flopy.utils.mflistfile` `ListBudget` class                   | **pandas** >= 0.15.0                                                     |
-| `.get_dataframes()` in `flopy.utils.observationfile` `ObsFiles` class                | **pandas** >= 0.15.0                                                     |
-| `.get_dataframes()` in `flopy.utils.sfroutputfile` `ModflowSfr2` class               | **pandas** >= 0.15.0                                                     |
-| `.get_dataframes()` in `flopy.utils.util_list` `MfList` class                        | **pandas** >= 0.15.0                                                     |
-| `.get_dataframes()` in `flopy.utils.zonebud` `ZoneBudget` class                      | **pandas** >= 0.15.0                                                     |
-| `.pivot_keyarray()` in `flopy.mf6.utils.arrayutils` `AdvancedPackageUtil` class      | **pandas** >= 0.15.0                                                     |
-| `._get_vertices()` in `flopy.mf6.utils.binaryfile_utils` `MFOutputRequester` class   | **pandas** >= 0.15.0                                                     |
-| `.get_dataframe()` in `flopy.mf6.utils.mfobservation` `Observations` class           | **pandas** >= 0.15.0                                                     |
-| `.df()` in `flopy.modflow.mfsfr2` `SfrFile` class                                    | **pandas** >= 0.15.0                                                     |
-| `.time_coverage()` in `flopy.export.metadata` `acc` class - ***used if available***  | **pandas** >= 0.15.0                                                     |
-| `.loadtxt()` in `flopy.utils.flopyio` - ***used if available***                      | **pandas** >= 0.15.0                                                     |
 | `.generate_classes()` in `flopy.mf6.utils`                                           | [**modflow-devtools**](https://github.com/MODFLOW-USGS/modflow-devtools) |
 | `GridIntersect()` in `flopy.utils.gridintersect`                                     | **shapely**                                                              |
 | `GridIntersect().plot_polygon()` in `flopy.utils.gridintersect`                      | **shapely** and **descartes**                                            |
diff --git a/flopy/export/metadata.py b/flopy/export/metadata.py
index d3325864fd..f7bc326857 100644
--- a/flopy/export/metadata.py
+++ b/flopy/export/metadata.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pandas as pd
 
 from ..utils import import_optional_dependency
 from ..utils.flopy_io import get_url_text
@@ -191,7 +192,6 @@ def time_coverage(self):
         -------
 
         """
-        pd = import_optional_dependency("pandas", errors="ignore")
 
         l = self.sb["dates"]
         tc = {}
diff --git a/flopy/mf6/utils/binaryfile_utils.py b/flopy/mf6/utils/binaryfile_utils.py
index 7ed85ed080..24fc9b8512 100644
--- a/flopy/mf6/utils/binaryfile_utils.py
+++ b/flopy/mf6/utils/binaryfile_utils.py
@@ -1,6 +1,7 @@
 import os
 
 import numpy as np
+import pandas as pd
 
 from ...utils import binaryfile as bf
 from ...utils import import_optional_dependency
@@ -231,11 +232,6 @@ def _get_vertices(mfdict, key):
         elevations corresponding to a row column location
         """
 
-        pd = import_optional_dependency(
-            "pandas",
-            error_message="MFOutputRequester._get_vertices() requires pandas.",
-        )
-
         mname = key[0]
         cellid = mfdict[(mname, "DISV8", "CELL2D", "cell2d_num")]
 
diff --git a/flopy/mf6/utils/mfobservation.py b/flopy/mf6/utils/mfobservation.py
index ab3fe30c3f..09a4c57081 100644
--- a/flopy/mf6/utils/mfobservation.py
+++ b/flopy/mf6/utils/mfobservation.py
@@ -1,6 +1,7 @@
 import csv
 
 import numpy as np
+import pandas as pd
 
 from ...utils import import_optional_dependency
 
@@ -210,10 +211,6 @@ def get_dataframe(
         pd.DataFrame
 
         """
-        pd = import_optional_dependency(
-            "pandas",
-            error_message="get_dataframe() requires pandas.",
-        )
 
         data_str = self._reader(self.Obsname)
         data = self._array_to_dict(data_str)
diff --git a/flopy/modflow/mfsfr2.py b/flopy/modflow/mfsfr2.py
index cf8355c01d..c2a1f43400 100644
--- a/flopy/modflow/mfsfr2.py
+++ b/flopy/modflow/mfsfr2.py
@@ -5,6 +5,7 @@
 import warnings
 
 import numpy as np
+import pandas as pd
 from numpy.lib import recfunctions
 
 from ..pakbase import Package
@@ -646,7 +647,6 @@ def paths(self):
 
     @property
     def df(self):
-        pd = import_optional_dependency("pandas")
         return pd.DataFrame(self.reach_data)
 
     def _make_graph(self):
@@ -1573,8 +1573,6 @@ def plot_path(self, start_seg=None, end_seg=0, plot_segment_lines=True):
         """
         import matplotlib.pyplot as plt
 
-        pd = import_optional_dependency("pandas")
-
         df = self.df
         m = self.parent
         mfunits = m.modelgrid.units
diff --git a/flopy/modpath/mp6sim.py b/flopy/modpath/mp6sim.py
index bf96a56cd1..1900b40c6d 100644
--- a/flopy/modpath/mp6sim.py
+++ b/flopy/modpath/mp6sim.py
@@ -8,6 +8,7 @@
 
 """
 import numpy as np
+import pandas as pd
 
 from ..pakbase import Package
 from ..utils import Util3d, import_optional_dependency
@@ -411,8 +412,8 @@ class StartingLocationsFile(Package):
         Input style described in MODPATH6 manual (currently only input style 1 is supported)
     extension : string
         Filename extension (default is 'loc')
-    use_pandas: bool, default False
-        If True and pandas is available use pandas to write the particle locations >2x speed
+    use_pandas: bool, default True
+        If True use pandas to write the particle locations >2x speed
     """
 
     def __init__(
@@ -421,7 +422,7 @@ def __init__(
         inputstyle=1,
         extension="loc",
         verbose=False,
-        use_pandas=False,
+        use_pandas=True,
     ):
         super().__init__(model, extension, "LOC", 33)
 
@@ -510,10 +511,6 @@ def _write_particle_data_with_pandas(self, data, float_format):
         :param save_group_mapper bool, if true, save a groupnumber to group name mapper as well.
         :return:
         """
-        pd = import_optional_dependency(
-            "pandas",
-            error_message="specify `use_pandas=False` to use slower methods without pandas",
-        )
         # convert float format string to pandas float format
         float_format = (
             float_format.replace("{", "").replace("}", "").replace(":", "%")
diff --git a/flopy/utils/flopy_io.py b/flopy/utils/flopy_io.py
index cfd3190d67..17078dabb7 100644
--- a/flopy/utils/flopy_io.py
+++ b/flopy/utils/flopy_io.py
@@ -9,6 +9,7 @@
 from typing import Union
 
 import numpy as np
+import pandas as pd
 
 
 def _fmt_string(array, float_format="{}"):
@@ -326,7 +327,7 @@ def loadtxt(
     file, delimiter=" ", dtype=None, skiprows=0, use_pandas=True, **kwargs
 ):
     """
-    Use pandas if it is available to load a text file
+    Use pandas to load a text file
     (significantly faster than n.loadtxt or genfromtxt see
     https://stackoverflow.com/q/18259393/)
 
@@ -352,15 +353,12 @@ def loadtxt(
     """
     from ..utils import import_optional_dependency
 
-    # test if pandas should be used, if available
     if use_pandas:
-        pd = import_optional_dependency("pandas")
         if delimiter.isspace():
             kwargs["delim_whitespace"] = True
         if isinstance(dtype, np.dtype) and "names" not in kwargs:
             kwargs["names"] = dtype.names
 
-    # if use_pandas and pd then use pandas
     if use_pandas:
         df = pd.read_csv(file, dtype=dtype, skiprows=skiprows, **kwargs)
         return df.to_records(index=False)
diff --git a/flopy/utils/mflistfile.py b/flopy/utils/mflistfile.py
index b88850d078..6ac362cef1 100644
--- a/flopy/utils/mflistfile.py
+++ b/flopy/utils/mflistfile.py
@@ -10,6 +10,7 @@
 import re
 
 import numpy as np
+import pandas as pd
 
 from ..utils import import_optional_dependency
 from ..utils.flopy_io import get_ts_sp
@@ -494,11 +495,6 @@ def get_dataframes(self, start_datetime="1-1-1970", diff=False):
 
         """
 
-        pd = import_optional_dependency(
-            "pandas",
-            error_message="ListBudget.get_dataframes() requires pandas.",
-        )
-
         if not self._isvalid:
             return None
         totim = self.get_times()
diff --git a/flopy/utils/mtlistfile.py b/flopy/utils/mtlistfile.py
index 5214859795..9f70c72b41 100644
--- a/flopy/utils/mtlistfile.py
+++ b/flopy/utils/mtlistfile.py
@@ -6,6 +6,7 @@
 import warnings
 
 import numpy as np
+import pandas as pd
 
 from ..utils import import_optional_dependency
 
@@ -76,11 +77,6 @@ def parse(
             (optionally) surface-water mass budget.
             If the SFT process is not used, df_sw is None.
         """
-        pd = import_optional_dependency(
-            "pandas",
-            error_message="MtListBudget.parse() requires pandas.",
-        )
-
         self.gw_data = {}
         self.sw_data = {}
         self.lcount = 0
@@ -182,11 +178,6 @@ def parse(
         return df_gw, df_sw
 
     def _diff(self, df):
-        pd = import_optional_dependency(
-            "pandas",
-            error_message="MtListBudget._diff() requires pandas.",
-        )
-
         out_cols = [
             c for c in df.columns if "_out" in c and not c.startswith("net_")
         ]
diff --git a/flopy/utils/observationfile.py b/flopy/utils/observationfile.py
index d385320daf..5f6f37fdfb 100644
--- a/flopy/utils/observationfile.py
+++ b/flopy/utils/observationfile.py
@@ -1,6 +1,7 @@
 import io
 
 import numpy as np
+import pandas as pd
 
 from ..utils import import_optional_dependency
 from ..utils.flopy_io import get_ts_sp
@@ -179,11 +180,6 @@ def get_dataframe(
 
         from ..utils.utils_def import totim_to_datetime
 
-        pd = import_optional_dependency(
-            "pandas",
-            error_message="ObsFiles.get_dataframe() requires pandas.",
-        )
-
         i0 = 0
         i1 = self.data.shape[0]
         if totim is not None:
diff --git a/flopy/utils/sfroutputfile.py b/flopy/utils/sfroutputfile.py
index 5e8d18ea72..797bcfdad6 100644
--- a/flopy/utils/sfroutputfile.py
+++ b/flopy/utils/sfroutputfile.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pandas as pd
 
 from ..utils import import_optional_dependency
 
@@ -52,8 +53,6 @@ def __init__(self, filename, geometries=None, verbose=False):
         Class constructor.
         """
 
-        self.pd = import_optional_dependency("pandas")
-
         # get the number of rows to skip at top, and the number of data columns
         self.filename = filename
         evaluated_format = False
@@ -172,13 +171,10 @@ def get_dataframe(self):
             "skiprows": self.sr,
             "low_memory": False,
         }
-        try:  # since pandas 1.3.0
-            df = self.pd.read_csv(**kwargs, on_bad_lines="skip")
-        except TypeError:  # before pandas 1.3.0
-            df = self.pd.read_csv(**kwargs, error_bad_lines=False)
+        df = pd.read_csv(**kwargs, on_bad_lines="skip")
 
         # drop text between stress periods; convert to numeric
-        df["layer"] = self.pd.to_numeric(df.layer, errors="coerce")
+        df["layer"] = pd.to_numeric(df.layer, errors="coerce")
         df.dropna(axis=0, inplace=True)
 
         # convert to proper dtypes
@@ -247,7 +243,7 @@ def get_results(self, segment, reach):
             results = self._get_result(segment, reach)
         except:
             locsr = list(zip(segment, reach))
-            results = self.pd.DataFrame()
+            results = pd.DataFrame()
             for s, r in locsr:
                 srresults = self._get_result(s, r)
                 if len(srresults) > 0:
diff --git a/flopy/utils/util_list.py b/flopy/utils/util_list.py
index 7b3b1b695a..2aa83be270 100644
--- a/flopy/utils/util_list.py
+++ b/flopy/utils/util_list.py
@@ -11,6 +11,7 @@
 import warnings
 
 import numpy as np
+import pandas as pd
 
 from ..datbase import DataInterface, DataListInterface, DataType
 from ..utils import import_optional_dependency
@@ -438,15 +439,7 @@ def get_dataframe(self, squeeze=False):
             stress periods where at least one cells is different,
             otherwise it is equal to the number of keys in MfList.data.
 
-        Notes
-        -----
-        Requires pandas.
-
         """
-        pd = import_optional_dependency(
-            "pandas",
-            error_message="MfList.get_dataframe() requires pandas.",
-        )
 
         # make a dataframe of all data for all stress periods
         names = ["per", "k", "i", "j"]
diff --git a/flopy/utils/utl_import.py b/flopy/utils/utl_import.py
index 321215c7a5..166102814c 100644
--- a/flopy/utils/utl_import.py
+++ b/flopy/utils/utl_import.py
@@ -49,7 +49,6 @@
 VERSIONS = {
     "shapefile": "2.0.0",
     "dateutil": "2.4.0",
-    "pandas": "0.15.0",
 }
 
 # A mapping from import name to package name (on PyPI) for packages where
diff --git a/flopy/utils/zonbud.py b/flopy/utils/zonbud.py
index 053b8e3df3..5740669c11 100644
--- a/flopy/utils/zonbud.py
+++ b/flopy/utils/zonbud.py
@@ -4,6 +4,7 @@
 from typing import Union
 
 import numpy as np
+import pandas as pd
 
 from . import import_optional_dependency
 from .utils_def import totim_to_datetime
@@ -2378,10 +2379,6 @@ def _recarray_to_dataframe(
 
     pd.DataFrame
     """
-    pd = import_optional_dependency(
-        "pandas",
-        error_message="ZoneBudget.get_dataframes() requires pandas.",
-    )
 
     valid_index_keys = ["totim", "kstpkper"]
     s = f'index_key "{index_key}" is not valid.'
@@ -3000,10 +2997,6 @@ def _volumetric_flux(recarray, modeltime, extrapolate_kper=False):
         pd.DataFrame
 
     """
-    pd = import_optional_dependency(
-        "pandas",
-        error_message="ZoneBudget._volumetric_flux() requires pandas.",
-    )
 
     nper = len(modeltime.nstp)
     volumetric_data = {}
diff --git a/pyproject.toml b/pyproject.toml
index 8eb30c8ef5..dac35baaed 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -30,6 +30,7 @@ requires-python = ">=3.8"
 dependencies = [
     "numpy >=1.15.0",
     "matplotlib >=1.4.0",
+    "pandas >=2.0.0"
 ]
 dynamic = ["version", "readme"]
 
@@ -65,7 +66,6 @@ optional = [
     "geojson",
     "imageio",
     "netcdf4",
-    "pandas",
     "pymetis ; platform_system != 'Windows'",
     "pyproj",
     "pyshp",