diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml
deleted file mode 100644
index 8f1af21105..0000000000
--- a/.github/workflows/benchmark.yml
+++ /dev/null
@@ -1,138 +0,0 @@
-name: Benchmarks
-
-on:
-  schedule:
-    - cron: '0 8 * * *' # run at 8 AM UTC (12 am PST)
-
-jobs:
-  benchmark:
-    name: Benchmarks
-    runs-on: ${{ matrix.os }}
-    strategy:
-      fail-fast: false
-      matrix:
-        os: [ ubuntu-latest, macos-latest, windows-latest ]
-        python-version: [ "3.10", "3.11", "3.12" ]
-    defaults:
-      run:
-        shell: bash -l {0}
-    timeout-minutes: 90
-
-    steps:
-      - name: Checkout repo
-        uses: actions/checkout@v6
-
-      - name: Setup Python ${{ matrix.python-version }}
-        uses: astral-sh/setup-uv@v7
-        with:
-          cache-dependency-glob: "**/pyproject.toml"
-          python-version: ${{ matrix.python-version }}
-
-      - name: Install FloPy
-        run: uv sync --all-extras
-
-      - name: Install Modflow executables
-        uses: modflowpy/install-modflow-action@v1
-
-      - name: Install triangle (macOS workaround)
-        if: runner.os == 'macOS'
-        uses: modflowpy/install-modflow-action@v1
-        with:
-          repo: executables
-          ostag: mac
-          subset: triangle
-
-      - name: Run benchmarks
-        working-directory: autotest
-        run: |
-          mkdir -p .benchmarks
-          uv run pytest -v --durations=0 --benchmark-only --benchmark-json .benchmarks/${{ matrix.os }}_python${{ matrix.python-version }}.json --keep-failed=.failed
-          ls .benchmarks
-        env:
-          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-
-      - name: Upload failed benchmark artifact
-        uses: actions/upload-artifact@v6
-        if: failure()
-        with:
-          name: failed-benchmark-${{ matrix.os }}-${{ matrix.python-version }}-${{ github.run_id }}
-          path: autotest/.failed/**
-
-      - name: Upload benchmark result artifact
-        uses: actions/upload-artifact@v6
-        with:
-          name: benchmarks-${{ matrix.os }}-${{ matrix.python-version }}-${{ github.run_id }}
-          path: autotest/.benchmarks/*.json
-          include-hidden-files: true
-
-  post_benchmark:
-    needs:
-      - benchmark
-    name: Process benchmark results
-    runs-on: ubuntu-latest
-    defaults:
-      run:
-        shell: bash
-    timeout-minutes: 10
-
-    steps:
-      - name: Checkout repo
-        uses: actions/checkout@v6
-
-      - name: Setup Python
-        uses: astral-sh/setup-uv@v7
-        with:
-          cache-dependency-glob: "**/pyproject.toml"
-
-      - name: Install FloPy
-        run: uv sync
-
-      - name: Install seaborn
-        run: uv pip install seaborn
-
-      - name: Download all artifacts
-        uses: actions/download-artifact@v7
-        with:
-          path: autotest/.benchmarks
-
-      - name: Process benchmark results
-        run: |
-          repo="${{ github.repository }}"
-          path="autotest/.benchmarks"
-
-          # list benchmark artifacts
-          artifact_json=$(gh api -X GET -H "Accept: application/vnd.github+json" /repos/$repo/actions/artifacts)
-
-          # get artifact ids and download artifacts
-          get_artifact_ids="
-          import json
-          import sys
-          from os import linesep
-
-          artifacts = json.load(sys.stdin, strict=False)['artifacts']
-          artifacts = [a for a in artifacts if a['name'].startswith('benchmarks-') and a['name'].split('-')[-1].isdigit()]
-
-          print(linesep.join([str(a['id']) for a in artifacts]))
-          "
-          echo $artifact_json \
-            | python -c "$get_artifact_ids" \
-            | xargs -I@ bash -c "gh api -H 'Accept: application/vnd.github+json' /repos/$repo/actions/artifacts/@/zip >> $path/@.zip"
-
-          # unzip artifacts
-          zipfiles=( $path/*.zip )
-          if (( ${#zipfiles[@]} )); then
-            unzip -o "$path/*.zip" -d $path
-          fi
-
-          # process benchmarks
-          uv run scripts/process_benchmarks.py $path $path
-        env:
-          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-
-      - name: Upload benchmark results
-        uses: actions/upload-artifact@v6
-        with:
-          name: benchmarks-${{ github.run_id }}
-          path: |
-            autotest/.benchmarks/*.csv
-            autotest/.benchmarks/*.png
diff --git a/.github/workflows/codspeed.yml b/.github/workflows/codspeed.yml
new file mode 100644
index 0000000000..2fc0b4cfe3
--- /dev/null
+++ b/.github/workflows/codspeed.yml
@@ -0,0 +1,78 @@
+name: CodSpeed Benchmarks
+
+on:
+  push:
+    branches:
+      - develop
+      - main
+  pull_request:
+  workflow_dispatch:
+
+# Required for OIDC authentication
+permissions:
+  contents: read
+  actions: read
+  id-token: write
+
+jobs:
+  benchmarks:
+    # Benchmarks are sharded by functionality:
+    # - input-io: model/simulation input file I/O
+    # - output-io: model output file readers
+    # - pre-post: pre/postprocessing, grids, rasters, arrays, export
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        shard:
+          - name: "input-io"
+            files: >-
+              benchmark_mf6_input.py
+              benchmark_mf2005_input.py
+          - name: "output-io"
+            files: >-
+              benchmark_cellbudgetfile.py
+              benchmark_zonebudget.py
+              benchmark_mf6listbudget.py
+              benchmark_headfile.py
+              benchmark_headufile.py
+              benchmark_formattedfile.py
+              benchmark_pathlinefile.py
+              benchmark_endpointfile.py
+              benchmark_mtlistfile.py
+              benchmark_sfroutputfile.py
+              benchmark_ucnfile.py
+              benchmark_mflistbudget.py
+              benchmark_mfusglistbudget.py
+          - name: "pre-post"
+            files: >-
+              benchmark_gridintersect.py
+              benchmark_grids.py
+              benchmark_rasters.py
+              benchmark_arrays.py
+              benchmark_export.py
+              benchmark_postprocessing.py
+    name: "benchmarks (${{ matrix.shard.name }})"
+    steps:
+      - name: Checkout repo
+        uses: actions/checkout@v4
+
+      - name: Setup Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+
+      - name: Install uv
+        uses: astral-sh/setup-uv@v4
+        with:
+          cache-dependency-glob: "**/pyproject.toml"
+
+      - name: Install FloPy with test dependencies
+        run: uv sync --all-extras
+
+      - name: Run benchmarks with CodSpeed
+        uses: CodSpeedHQ/action@v3
+        with:
+          run: cd autotest/benchmarks && uv run pytest ${{ matrix.shard.files }} --codspeed
+        env:
+          CODSPEED_SHARD_NAME: ${{ matrix.shard.name }}
diff --git a/DEVELOPER.md b/DEVELOPER.md
index 0f0943a75e..f7a456bace 100644
--- a/DEVELOPER.md
+++ b/DEVELOPER.md
@@ -370,45 +370,81 @@ To allow optional separation of performance from correctness concerns, performan
 
 #### Benchmarking
 
-Any test function can be turned into a benchmark by requesting the `benchmark` fixture (i.e. declaring a `benchmark` argument), which can be used to wrap any function call. For instance:
+FloPy includes a benchmark suite to track performance over time and validate optimization efforts. Benchmarking focuses on I/O operations, data structure manipulations, and utility functions.
 
-```python
-def test_benchmark(benchmark):
-    def sleep_1s():
-        import time
-        time.sleep(1)
-        return True
+**Note**: Benchmarks test FloPy code performance only, not the runtime of the various executables FloPy drives.
+
+Benchmarks use [CodSpeed](https://codspeed.io) to automatically track performance in CI. Benchmarks written using `pytest-benchmark` syntax are compatible. Benchmarks are organized in `autotest/benchmarks/` by functional area.
+
+##### Running Benchmarks
+
+```bash
+# Run all benchmarks
+pytest autotest/benchmarks --benchmark-only
+
+# Run specific benchmark file
+pytest autotest/benchmarks/benchmark_io_mf6.py --benchmark-only
 
-    assert benchmark(sleep_1s)
+# Run benchmarks by marker
+pytest -m "benchmark and not slow" --benchmark-only
+
+# Save results to file
+pytest autotest/benchmarks --benchmark-only --benchmark-autosave
+
+# Compare against saved baseline
+pytest autotest/benchmarks --benchmark-only --benchmark-compare
 ```
 
-Arguments can be provided to the function as well:
+##### Writing Benchmarks
+
+Any test function can be turned into a benchmark by requesting the `benchmark` fixture:
 
 ```python
-def test_benchmark(benchmark):
-    def sleep_s(s):
-        import time
-        time.sleep(s)
-        return True
+@pytest.mark.benchmark
+def test_model_load_time(benchmark, function_tmpdir):
+    """
+    Benchmark model loading time.
+
+    Measures time to load a MODFLOW model from disk.
+    """
+    model = create_test_model(function_tmpdir)
+    model.write_input()
 
-    assert benchmark(sleep_s, 1)
+    benchmark(lambda: Modflow.load(f"{model.name}.nam", model_ws=function_tmpdir))
 ```
 
-Rather than alter an existing function call to use this syntax, a lambda can be used to wrap the call unmodified:
+**Best Practices:**
+- Use descriptive test names (e.g., `test_mf6_sim_load_large`, not `test_load1`)
+- Include docstrings explaining what is benchmarked and why
+- Use fixtures for setup (not timed)
+- Mark all benchmarks `@pytest.mark.benchmark`
+- Mark slow benchmarks with `@pytest.mark.slow`
+
+##### Advanced Usage
+
+Arguments can be provided to benchmark functions:
 
 ```python
-def test_benchmark(benchmark):
-    def sleep_s(s):
-        import time
-        time.sleep(s)
-        return True
+def test_benchmark_with_args(benchmark):
+    benchmark(some_function, arg1, arg2)
+```
+
+For fine-grained control over iterations and rounds:
 
-    assert benchmark(lambda: sleep_s(1))
+```python
+def test_benchmark_controlled(benchmark):
+    benchmark.pedantic(some_function, iterations=10, rounds=5)
 ```
 
-This can be convenient when the function call is complicated or passes many arguments.
+Lambda functions are convenient for wrapping complex calls:
+
+```python
+def test_complex_benchmark(benchmark):
+    result = benchmark(lambda: complex_function(many, different, args))
+    assert result is not None
+```
 
-Benchmarked functions are repeated several times (the number of iterations depending on the test's runtime, with faster tests generally getting more reps) to compute summary statistics. To control the number of repetitions and rounds (repetitions of repetitions) use `benchmark.pedantic`, e.g. `benchmark.pedantic(some_function(), iterations=1, rounds=1)`.
+##### Configuration
 
 Benchmarking is incompatible with `pytest-xdist` and is disabled automatically when tests are run in parallel. When tests are not run in parallel, benchmarking is enabled by default. Benchmarks can be disabled with the `--benchmark-disable` flag.
 
diff --git a/autotest/benchmarks/__init__.py b/autotest/benchmarks/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/autotest/benchmarks/benchmark_arrays.py b/autotest/benchmarks/benchmark_arrays.py
new file mode 100644
index 0000000000..8678093676
--- /dev/null
+++ b/autotest/benchmarks/benchmark_arrays.py
@@ -0,0 +1,104 @@
+"""
+Benchmarks for flopy.utils.Util2d and Util3d operations including:
+- Array creation
+- Binary file I/O
+- Array access performance
+"""
+
+import numpy as np
+import pytest
+
+from flopy.modflow import Modflow
+from flopy.utils import Util2d, Util3d
+
+SIZES = {
+    "small": {"nlay": 3, "nrow": 10, "ncol": 10},
+    "medium": {"nlay": 10, "nrow": 100, "ncol": 100},
+    "large": {"nlay": 20, "nrow": 200, "ncol": 200},
+}
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("size", ["small", "medium", "large"])
+def test_util2d_create(benchmark, function_tmpdir, size):
+    dims = SIZES[size]
+    shape = (dims["nrow"], dims["ncol"])
+    data = np.random.random(shape)
+    ml = Modflow(model_ws=function_tmpdir)
+
+    def create_util2d():
+        return Util2d(ml, shape, np.float32, data.copy(), "test")
+
+    benchmark(create_util2d)
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("size", ["small", "medium", "large"])
+def test_util3d_create(benchmark, function_tmpdir, size):
+    dims = SIZES[size]
+    shape = (dims["nlay"], dims["nrow"], dims["ncol"])
+    data = np.random.random(shape)
+    ml = Modflow(model_ws=function_tmpdir)
+
+    def create_util3d():
+        return Util3d(ml, shape, np.float32, data.copy(), "test")
+
+    benchmark(create_util3d)
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("size", ["small", "medium", "large"])
+def test_util2d_external_write(benchmark, function_tmpdir, size):
+    dims = SIZES[size]
+    shape = (dims["nrow"], dims["ncol"])
+    data = np.random.random(shape).astype(np.float32)
+    fpath = function_tmpdir / "test_array.dat"
+
+    def write_bin():
+        Util2d.write_bin(shape, fpath, data, bintype="head")
+
+    benchmark(write_bin)
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("size", ["small", "medium", "large"])
+def test_util3d_external_write(benchmark, function_tmpdir, size):
+    dims = SIZES[size]
+    shape = (dims["nlay"], dims["nrow"], dims["ncol"])
+    data = np.random.random(shape).astype(np.float32)
+    fpath = function_tmpdir / "test_array3d.dat"
+
+    def write_bin():
+        for i in range(shape[0]):
+            layer_path = function_tmpdir / f"test_array3d_lay{i}.dat"
+            Util2d.write_bin((shape[1], shape[2]), layer_path, data[i], bintype="head")
+
+    benchmark(write_bin)
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("size", ["small", "medium", "large"])
+def test_util2d_array_copy(benchmark, function_tmpdir, size):
+    dims = SIZES[size]
+    shape = (dims["nrow"], dims["ncol"])
+    data = np.random.random(shape)
+    ml = Modflow(model_ws=function_tmpdir)
+    u2d = Util2d(ml, shape, np.float32, data, "test")
+    benchmark(lambda: u2d.array.copy())
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("size", ["small", "medium", "large"])
+def test_util3d_array_copy(benchmark, function_tmpdir, size):
+    dims = SIZES[size]
+    shape = (dims["nlay"], dims["nrow"], dims["ncol"])
+    data = np.random.random(shape)
+    ml = Modflow(model_ws=function_tmpdir)
+    u3d = Util3d(ml, shape, np.float32, data, "test")
+    benchmark(lambda: u3d.array.copy())
diff --git a/autotest/benchmarks/benchmark_cellbudgetfile.py b/autotest/benchmarks/benchmark_cellbudgetfile.py
new file mode 100644
index 0000000000..1ae3226347
--- /dev/null
+++ b/autotest/benchmarks/benchmark_cellbudgetfile.py
@@ -0,0 +1,54 @@
+from pathlib import Path
+
+import pytest
+
+from flopy.utils import CellBudgetFile
+
+
+@pytest.fixture
+def cbcf(example_data_path) -> CellBudgetFile:
+    return CellBudgetFile(
+        example_data_path
+        / "mf6"
+        / "create_tests"
+        / "test021_twri"
+        / "expected_output"
+        / "twri.cbc"
+    )
+
+
+@pytest.mark.benchmark
+def test_cellbudgetfile_load(benchmark, cbcf):
+    benchmark(lambda: CellBudgetFile(cbcf.filename))
+
+
+@pytest.mark.benchmark
+def test_cellbudgetfile_get_data_all(benchmark, cbcf):
+    unique_records = cbcf.headers[["text", "imeth"]].drop_duplicates()
+    term = unique_records.iloc[0]["text"]
+    benchmark(lambda: cbcf.get_data(text=term))
+
+
+@pytest.mark.benchmark
+def test_cellbudgetfile_get_data_one(benchmark, cbcf):
+    benchmark(lambda: cbcf.get_data(kstpkper=(0, 1)))
+
+
+@pytest.mark.benchmark
+def test_cellbudgetfile_list_records(benchmark, cbcf):
+    benchmark(cbcf.list_records)
+
+
+@pytest.mark.benchmark
+def test_cellbudgetfile_list_unique_records(benchmark, cbcf):
+    benchmark(cbcf.list_unique_records)
+
+
+@pytest.mark.benchmark
+def test_cellbudgetfile_get_times(benchmark, cbcf):
+    benchmark(cbcf.get_times)
+
+
+@pytest.mark.benchmark
+def test_cellbudgetfile_get_kstpkper(benchmark, cbcf):
+    benchmark(cbcf.get_kstpkper)
diff --git a/autotest/benchmarks/benchmark_endpointfile.py b/autotest/benchmarks/benchmark_endpointfile.py
new file mode 100644
index 0000000000..1f2b8f4100
--- /dev/null
+++ b/autotest/benchmarks/benchmark_endpointfile.py
@@ -0,0 +1,37 @@
+import pytest
+
+from autotest.benchmarks.benchmark_pathlinefile import ex01_mp7_model
+from autotest.test_mp7 import ex01_mf6_model
+from flopy.utils.modpathfile import EndpointFile
+
+
+@pytest.fixture
+def epf(ex01_mp7_model) -> EndpointFile:
+    mp, ws = ex01_mp7_model
+    mp.write_input()
+    success, buff = mp.run_model()
+    assert success
+    return EndpointFile(ws / f"{mp.name}.mpend")
+
+
+@pytest.mark.benchmark
+def test_endpointfile_load(benchmark, epf):
+    benchmark(lambda: EndpointFile(epf.fname))
+
+
+@pytest.mark.benchmark
+def test_endpointfile_get_data(benchmark, epf):
+    benchmark(epf.get_data)
+
+
+@pytest.mark.benchmark
+def test_endpointfile_get_alldata(benchmark, epf):
+    benchmark(epf.get_alldata)
+
+
+@pytest.mark.benchmark
+def test_endpointfile_to_geodataframe(benchmark, ex01_mf6_model, epf):
+    pytest.importorskip("geopandas")
+    sim, function_tmpdir = ex01_mf6_model
+    gwf = sim.get_model()
+    benchmark(lambda: epf.to_geodataframe(gwf.modelgrid))
diff --git a/autotest/benchmarks/benchmark_export.py b/autotest/benchmarks/benchmark_export.py
new file mode 100644
index 0000000000..bcdbe44bfe
--- /dev/null
+++ b/autotest/benchmarks/benchmark_export.py
@@ -0,0 +1,52 @@
+import pytest
+from modflow_devtools.misc import has_pkg
+
+from .conftest import load_mf6_sim, load_mf2005_model
+
+
+@pytest.mark.benchmark
+def test_mf6_export_shapefile(benchmark, function_tmpdir):
+    sim = load_mf6_sim(function_tmpdir, model_key="freyberg")
+    gwf = sim.get_model()
+    output_path = function_tmpdir / "export.shp"
+    benchmark(lambda: gwf.export(str(output_path)))
+
+
+@pytest.mark.benchmark
+def test_mf2005_export_shapefile(benchmark, function_tmpdir):
+    model = load_mf2005_model(function_tmpdir, model_key="freyberg")
+    output_path = function_tmpdir / "export_mf2005.shp"
+    benchmark(lambda: model.export(str(output_path)))
+
+
+@pytest.mark.benchmark
+@pytest.mark.skipif(not has_pkg("netCDF4"), reason="requires netCDF4")
+def test_mf6_export_netcdf(benchmark, function_tmpdir):
+    import uuid
+
+    sim = load_mf6_sim(function_tmpdir, model_key="freyberg")
+    gwf = sim.get_model()
+
+    def export_netcdf():
+        # Use unique filename for each iteration to avoid file locking issues on Windows
+        output_path = function_tmpdir / f"export_{uuid.uuid4().hex[:8]}.nc"
+        gwf.export(str(output_path), fmt="netcdf")
+
+    benchmark(export_netcdf)
+
+
+@pytest.mark.benchmark
+@pytest.mark.skipif(not has_pkg("geopandas"), reason="requires geopandas")
+def test_mf6_modelgrid_to_geodataframe(benchmark, function_tmpdir):
+    sim = load_mf6_sim(function_tmpdir, model_key="freyberg")
+    gwf = sim.get_model()
+    benchmark(gwf.modelgrid.to_geodataframe)
+
+
+@pytest.mark.benchmark
+@pytest.mark.skipif(not has_pkg("vtk"), reason="requires vtk")
+def test_mf6_export_vtk(benchmark, function_tmpdir):
+    sim = load_mf6_sim(function_tmpdir, model_key="freyberg")
+    gwf = sim.get_model()
+    output_path = function_tmpdir / "export.vtk"
+    benchmark(lambda: gwf.export(str(output_path), fmt="vtk"))
diff --git a/autotest/benchmarks/benchmark_formattedfile.py b/autotest/benchmarks/benchmark_formattedfile.py
new file mode 100644
index 0000000000..b57aff8d8e
--- /dev/null
+++ b/autotest/benchmarks/benchmark_formattedfile.py
@@ -0,0 +1,70 @@
+"""
+Benchmarks for flopy.utils.formattedfile.FormattedHeadFile operations including:
+- Formatted (ASCII) head file loading
+- Data extraction for single stress period
+- Time series extraction
+- Full data retrieval
+"""
+
+import numpy as np
+import pytest
+
+from flopy.utils.formattedfile import FormattedHeadFile
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+def test_formattedfile_load(benchmark, example_data_path):
+    pth = example_data_path / "mf2005_test" / "test1tr.githds"
+    benchmark(lambda: FormattedHeadFile(pth))
+
+
+@pytest.fixture
+def fhd(example_data_path) -> FormattedHeadFile:
+    return FormattedHeadFile(str(example_data_path / "mf2005_test" / "test1tr.githds"))
+
+
+@pytest.mark.benchmark
+def test_formattedfile_get_data_totim(benchmark, fhd):
+    times = fhd.get_times()
+    mid_time = times[len(times) // 2]
+    benchmark(lambda: fhd.get_data(totim=mid_time))
+
+
+@pytest.mark.benchmark
+def test_formattedfile_get_data_first(benchmark, fhd):
+    benchmark(lambda: fhd.get_data(idx=0))
+
+
+@pytest.mark.benchmark
+def test_formattedfile_get_data_last(benchmark, fhd):
+    times = fhd.get_times()
+    benchmark(lambda: fhd.get_data(totim=times[-1]))
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+def test_formattedfile_get_alldata(benchmark, fhd):
+    benchmark(fhd.get_alldata)
+
+
+@pytest.mark.benchmark
+def test_formattedfile_get_ts(benchmark, fhd):
+    # Use a valid cell index based on test file dimensions (1, 15, 10)
+    benchmark(lambda: fhd.get_ts((0, 7, 5)))
+
+
+@pytest.mark.benchmark
+def test_formattedfile_get_times(benchmark, fhd):
+    benchmark(fhd.get_times)
+
+
+@pytest.mark.benchmark
+def test_formattedfile_get_kstpkper(benchmark, fhd):
+    benchmark(fhd.get_kstpkper)
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+def test_formattedfile_many_stress_periods(benchmark, fhd):
+    benchmark(fhd.get_alldata)
diff --git a/autotest/benchmarks/benchmark_gridintersect.py b/autotest/benchmarks/benchmark_gridintersect.py
new file mode 100644
index 0000000000..b4de7a902a
--- /dev/null
+++ b/autotest/benchmarks/benchmark_gridintersect.py
@@ -0,0 +1,183 @@
+"""
+Benchmarks for spatial intersection operations including:
+
+1. GridIntersect class (flopy.utils.gridintersect.GridIntersect):
+   - Initialization (with/without STR-tree spatial index)
+   - Point, LineString, and Polygon intersections
+   - Spatial query performance
+
+2. Grid.intersect() method (direct coordinate-based intersection):
+   - Single point lookup
+   - Batch point operations
+   - 3D coordinate intersection
+"""
+
+import numpy as np
+import pytest
+
+from autotest.test_grid_cases import GridCases
+from flopy.utils.geometry import LineString, Point, Polygon
+from flopy.utils.gridintersect import GridIntersect
+
+STRUCTURED_GRIDS = {
+    "small": GridCases.structured_small(),
+    "medium": GridCases.structured_medium(),
+    "large": GridCases.structured_large(),
+}
+
+
+# GridIntersect class benchmarks
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+@pytest.mark.parametrize("rtree", [True, False], ids=["rtree", "no_rtree"])
+def test_init(benchmark, grid, rtree):
+    benchmark(lambda: GridIntersect(grid, rtree=rtree))
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+@pytest.mark.parametrize("rtree", [True, False], ids=["rtree", "no_rtree"])
+def test_intersect_point(benchmark, grid, rtree):
+    gi = GridIntersect(grid, rtree=rtree)
+    xmin, xmax, ymin, ymax = grid.extent
+    point = Point((xmin + xmax) / 2, (ymin + ymax) / 2)
+    benchmark(lambda: gi.intersect(point, "point"))
+
+
+def make_line(grid, line_type) -> LineString:
+    xmin, xmax, ymin, ymax = grid.extent
+    if line_type == "diagonal":
+        return LineString([(xmin, ymin), (xmax, ymax)])
+    elif line_type == "horizontal":
+        y_mid = (ymin + ymax) / 2
+        return LineString([(xmin, y_mid), (xmax, y_mid)])
+    elif line_type == "complex":
+        x = np.linspace(xmin, xmax, 20)
+        y_mid = (ymin + ymax) / 2
+        y_range = (ymax - ymin) * 0.2
+        y = y_mid + y_range * np.sin(x / (xmax - xmin) * 10)
+        coords = list(zip(x, y))
+        return LineString(coords)
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+@pytest.mark.parametrize("line", ["diagonal", "horizontal", "complex"])
+@pytest.mark.parametrize("rtree", [True, False], ids=["rtree", "no_rtree"])
+def test_intersect_linestring_diagonal(benchmark, grid, line, rtree):
+    gi = GridIntersect(grid, rtree=rtree)
+    line = make_line(grid, line)
+    benchmark(lambda: gi.intersect(line, "linestring"))
+
+
+def make_poly(grid, poly_type) -> Polygon:
+    xmin, xmax, ymin, ymax = grid.extent
+    x_center = (xmin + xmax) / 2
+    y_center = (ymin + ymax) / 2
+    x_range = xmax - xmin
+    y_range = ymax - ymin
+
+    if poly_type == "small":
+        # 10% of grid extent centered
+        dx = x_range * 0.05
+        dy = y_range * 0.05
+        coords = [
+            (x_center - dx, y_center - dy),
+            (x_center + dx, y_center - dy),
+            (x_center + dx, y_center + dy),
+            (x_center - dx, y_center + dy),
+        ]
+    elif poly_type == "medium":
+        # 50% of grid extent centered
+        dx = x_range * 0.25
+        dy = y_range * 0.25
+        coords = [
+            (x_center - dx, y_center - dy),
+            (x_center + dx, y_center - dy),
+            (x_center + dx, y_center + dy),
+            (x_center - dx, y_center + dy),
+        ]
+    elif poly_type == "large":
+        # 90% of grid extent centered
+        dx = x_range * 0.45
+        dy = y_range * 0.45
+        coords = [
+            (x_center - dx, y_center - dy),
+            (x_center + dx, y_center - dy),
+            (x_center + dx, y_center + dy),
+            (x_center - dx, y_center + dy),
+        ]
+    elif poly_type == "irregular":
+        # Irregular star-like polygon
+        n_points = 20
+        theta = np.linspace(0, 2 * np.pi, n_points, endpoint=False)
+        r_base = min(x_range, y_range) * 0.3
+        r_var = r_base * 0.33
+        r = r_base + r_var * np.sin(5 * theta)
+        x = x_center + r * np.cos(theta)
+        y = y_center + r * np.sin(theta)
+        coords = list(zip(x, y))
+
+    return Polygon(coords)
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+@pytest.mark.parametrize("poly", ["small", "medium", "large", "irregular"])
+@pytest.mark.parametrize("rtree", [True, False], ids=["rtree", "no_rtree"])
+def test_intersect_polygon(benchmark, grid, poly, rtree):
+    gi = GridIntersect(grid, rtree=rtree)
+    polygon = make_poly(grid, poly)
+    benchmark(lambda: gi.intersect(polygon, "polygon"))
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+@pytest.mark.parametrize("poly", ["small", "medium", "large", "irregular"])
+@pytest.mark.parametrize("rtree", [True, False], ids=["rtree", "no_rtree"])
+def test_query_grid_polygon(benchmark, grid, poly, rtree):
+    gi = GridIntersect(grid, rtree=rtree)
+    polygon = make_poly(grid, poly)
+    benchmark(lambda: gi.query_grid(polygon))
+
+
+# Grid.intersect() method benchmarks (coordinate-based intersection)
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+def test_grid_intersect_single_point(benchmark, grid):
+    xmin, xmax, ymin, ymax = grid.extent
+    x_center = (xmin + xmax) / 2
+    y_center = (ymin + ymax) / 2
+    benchmark(lambda: grid.intersect(x_center, y_center))
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+def test_grid_intersect_batch_points(benchmark, grid):
+    xmin, xmax, ymin, ymax = grid.extent
+    x = np.linspace(xmin + 0.1 * (xmax - xmin), xmax - 0.1 * (xmax - xmin), 100)
+    y = np.linspace(ymin + 0.1 * (ymax - ymin), ymax - 0.1 * (ymax - ymin), 100)
+    xx, yy = np.meshgrid(x, y)
+    benchmark(lambda: grid.intersect(xx.ravel(), yy.ravel()))
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+def test_grid_intersect_3d(benchmark, grid):
+    xmin, xmax, ymin, ymax = grid.extent
+    x = np.linspace(xmin + 0.1 * (xmax - xmin), xmax - 0.1 * (xmax - xmin), 50)
+    y = np.linspace(ymin + 0.1 * (ymax - ymin), ymax - 0.1 * (ymax - ymin), 50)
+    xx, yy = np.meshgrid(x, y)
+    zz = np.ones_like(xx) * 5.0
+    benchmark(lambda: grid.intersect(xx.ravel(), yy.ravel(), zz.ravel()))
diff --git a/autotest/benchmarks/benchmark_grids.py b/autotest/benchmarks/benchmark_grids.py
new file mode 100644
index 0000000000..924ec1474f
--- /dev/null
+++ b/autotest/benchmarks/benchmark_grids.py
@@ -0,0 +1,67 @@
+"""
+Benchmarks for flopy.discretization grid operations including:
+- cellid <-> node number conversions
+- grid intersection operations
+- grid geometry properties
+"""
+
+import numpy as np
+import pytest
+
+from autotest.test_grid_cases import GridCases
+
+STRUCTURED_GRIDS = {
+    "small": GridCases.structured_small(),
+    "medium": GridCases.structured_medium(),
+    "large": GridCases.structured_large(),
+}
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+def test_structured_grid_get_lrc(benchmark, grid):
+    nodes = list(range(grid.nnodes))
+    benchmark(lambda: grid.get_lrc(nodes))
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+def test_structured_grid_get_node(benchmark, grid):
+    cellids = [
+        (lay, row, col)
+        for lay in range(grid.nlay)
+        for row in range(0, grid.nrow)
+        for col in range(0, grid.ncol)
+    ]
+    benchmark(lambda: grid.get_node(cellids))
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+def test_structured_grid_intersect_linestring(benchmark, grid):
+    # Create x, y coordinates along a diagonal line across the grid
+    x = np.linspace(0, grid.ncol, 100)
+    y = np.linspace(0, grid.nrow, 100)
+    benchmark(lambda: grid.intersect(x, y, forgive=True))
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+def test_structured_grid_intersect_point(benchmark, grid):
+    # Use grid center point
+    x = grid.ncol / 2.0
+    y = grid.nrow / 2.0
+    benchmark(lambda: grid.intersect(x, y))
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+def test_structured_grid_extent(benchmark, grid):
+    benchmark(lambda: grid.extent)
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize("grid", STRUCTURED_GRIDS.values(), ids=STRUCTURED_GRIDS.keys())
+def test_structured_grid_xyzcellcenters(benchmark, grid):
+    benchmark(lambda: grid.xyzcellcenters)
diff --git a/autotest/benchmarks/benchmark_headfile.py b/autotest/benchmarks/benchmark_headfile.py
new file mode 100644
index 0000000000..c38805a3d5
--- /dev/null
+++ b/autotest/benchmarks/benchmark_headfile.py
@@ -0,0 +1,83 @@
+"""
+Benchmarks for flopy.utils.HeadFile operations including:
+- HeadFile initialization
+- get_data() for single time steps
+- get_alldata() for full file reading
+- get_ts() for time series extraction
+"""
+
+from pathlib import Path
+
+import pytest
+
+from flopy.utils import HeadFile
+
+
+@pytest.mark.benchmark
+def test_headfile_load(benchmark, example_data_path):
+    pth = (
+        example_data_path
+        / "mf6"
+        / "create_tests"
+        / "test021_twri"
+        / "expected_output"
+        / "twri.hds"
+    )
+    benchmark(lambda: HeadFile(pth))
+
+
+@pytest.fixture
+def hdsf(example_data_path) -> HeadFile:
+    return HeadFile(
+        example_data_path
+        / "mf6"
+        / "create_tests"
+        / "test021_twri"
+        / "expected_output"
+        / "twri.hds"
+    )
+
+
+@pytest.mark.benchmark
+def test_headfile_get_data_single(benchmark, hdsf):
+    times = hdsf.get_times()
+    mid_time = times[len(times) // 2] if len(times) > 1 else times[0]
+    benchmark(lambda: hdsf.get_data(totim=mid_time))
+
+
+@pytest.mark.benchmark
+def test_headfile_get_alldata(benchmark, hdsf):
+    benchmark(hdsf.get_alldata)
+
+
+@pytest.mark.benchmark
+def test_headfile_get_ts(benchmark, hdsf):
+    benchmark(lambda: hdsf.get_ts((0, 10, 10)))
+
+
+@pytest.mark.benchmark
+def test_headfile_get_kstpkper(benchmark, hdsf):
+    # Use the first available kstpkper from the file
+    kstpkpers = hdsf.get_kstpkper()
+    kstpkper = kstpkpers[0] if kstpkpers else (0, 0)
+    benchmark(lambda: hdsf.get_data(kstpkper=kstpkper))
+
+
+@pytest.mark.benchmark
+def test_headfile_list_records(benchmark, hdsf):
+    benchmark(hdsf.list_records)
+
+
+@pytest.mark.benchmark
+def test_headfile_get_times(benchmark, hdsf):
+    benchmark(hdsf.get_times)
+
+
+@pytest.mark.benchmark
+def test_headfile_get_kstpkper_list(benchmark, hdsf):
+    benchmark(hdsf.get_kstpkper)
+
+
+@pytest.mark.benchmark
+def test_headfile_get_alldata_mf6(benchmark, hdsf):
+    benchmark(hdsf.get_alldata)
diff --git a/autotest/benchmarks/benchmark_headufile.py b/autotest/benchmarks/benchmark_headufile.py
new file mode 100644
index 0000000000..f1efa814f7
--- /dev/null
+++ b/autotest/benchmarks/benchmark_headufile.py
@@ -0,0 +1,32 @@
+import pytest
+
+from flopy.utils import HeadUFile
+
+
+@pytest.mark.benchmark
+def test_headufile_load(benchmark, example_data_path):
+    hds_file = example_data_path / "unstructured" / "headu.githds"
+    benchmark(lambda: HeadUFile(str(hds_file)))
+
+
+@pytest.fixture
+def huf(example_data_path) -> HeadUFile:
+    return HeadUFile(str(example_data_path / "unstructured" / "headu.githds"))
+
+
+@pytest.mark.benchmark
+def test_headufile_get_data(benchmark, huf):
+    times = huf.get_times()
+    mid_time = times[len(times) // 2] if len(times) > 0 else times[0]
+    benchmark(lambda: huf.get_data(totim=mid_time))
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+def test_headufile_get_alldata(benchmark, huf):
+    benchmark(huf.get_alldata)
+
+
+@pytest.mark.benchmark
+def test_headufile_get_ts(benchmark, huf):
+    benchmark(lambda: huf.get_ts(0))
diff --git a/autotest/benchmarks/benchmark_mf2005_input.py b/autotest/benchmarks/benchmark_mf2005_input.py
new file mode 100644
index 0000000000..95a11c81d2
--- /dev/null
+++ b/autotest/benchmarks/benchmark_mf2005_input.py
@@ -0,0 +1,62 @@
+from shutil import copytree
+
+import pytest
+from modflow_devtools.misc import get_namefile_paths
+
+from flopy.modflow import Modflow
+
+from .conftest import get_examples_path
+
+# Model directories from examples/data
+MODEL_DIRS = [
+    "freyberg_multilayer_transient",
+    "mf2005_test",
+]
+
+
+def pytest_generate_tests(metafunc):
+    """Dynamically parametrize tests with model names."""
+    if "model_name" in metafunc.fixturenames:
+        examples_path = get_examples_path()
+        models = []
+
+        # Only include models that exist
+        for model_dir in MODEL_DIRS:
+            if (examples_path / model_dir).exists():
+                models.append(model_dir)
+
+        metafunc.parametrize("model_name", models, ids=models)
+
+
+def _load_model(ws, model_name):
+    """Load MODFLOW-2005 model from examples."""
+    examples_path = get_examples_path()
+    source_path = examples_path / model_name
+    copytree(source_path, ws, dirs_exist_ok=True)
+
+    nam_files = get_namefile_paths(source_path, namefile="*.nam")
+    assert nam_files, f"No .nam file found in {source_path}"
+    nam_file = nam_files[0].name
+
+    return Modflow.load(nam_file, model_ws=ws, check=False)
+
+
+@pytest.mark.benchmark
+def test_mf2005_load(benchmark, function_tmpdir, model_name):
+    benchmark(lambda: _load_model(function_tmpdir, model_name))
+
+
+@pytest.mark.benchmark
+def test_mf2005_write_freyberg(benchmark, function_tmpdir):
+    ml = _load_model(function_tmpdir, "freyberg_multilayer_transient")
+    benchmark(ml.write_input)
+
+
+@pytest.mark.benchmark
+def test_mf2005_round_trip_freyberg(benchmark, function_tmpdir):
+    def round_trip():
+        ml = _load_model(function_tmpdir, "freyberg_multilayer_transient")
+        ml.write_input()
+        return Modflow.load("freyberg.nam", model_ws=function_tmpdir, check=False)
+
+    benchmark(round_trip)
diff --git a/autotest/benchmarks/benchmark_mf6_input.py b/autotest/benchmarks/benchmark_mf6_input.py
new file mode 100644
index 0000000000..2380470e4c
--- /dev/null
+++ b/autotest/benchmarks/benchmark_mf6_input.py
@@ -0,0 +1,77 @@
+"""
+Benchmarks for MODFLOW 6 I/O operations using models from:
+- examples/data (local models)
+- modflow-devtools registry (downloaded on-demand)
+"""
+
+from pathlib import Path
+from shutil import copytree
+
+import pytest
+from modflow_devtools.models import DEFAULT_REGISTRY, LocalRegistry
+
+from flopy.mf6 import MFSimulation
+
+from .conftest import get_examples_path
+
+# prefixes into the model registry
+PREFIXES = ["mf6/test", "mf6/large", "mf2005"]
+
+
+def pytest_generate_tests(metafunc):
+    # Use the --models-path command line option once or more to specify
+    # model directories. If at least one --models_path is provided,
+    # external tests (i.e. those using models from an external repo)
+    # will run against model input files found in the given location
+    # on the local filesystem rather than model input files from the
+    # official model registry. This is useful for testing changes to
+    # test model input files during MF6 development. See conftest.py
+    # for the models_path fixture and CLI argument definitions.
+    if "model_name" in metafunc.fixturenames:
+        # Try to get the models-path option, default to None if not available
+        try:
+            models_paths = metafunc.config.getoption("--models-path") or []
+        except ValueError:
+            models_paths = []
+
+        models_paths = [Path(p).expanduser().resolve().absolute() for p in models_paths]
+        registry = LocalRegistry() if any(models_paths) else DEFAULT_REGISTRY
+        registry_type = type(registry).__name__.lower().replace("registry", "")
+        metafunc.parametrize("registry", [registry], ids=[registry_type])
+        models = []
+        if "local" in registry_type:
+            try:
+                namefile_pattern = (
+                    metafunc.config.getoption("--namefile-pattern") or "mfsim.nam"
+                )
+            except ValueError:
+                namefile_pattern = "mfsim.nam"
+            for path in models_paths:
+                registry.index(path, namefile=namefile_pattern)
+            models.extend(registry.models.keys())
+        else:
+            for model_prefix in PREFIXES:
+                models.extend(
+                    [m for m in registry.models.keys() if m.startswith(model_prefix)]
+                )
+        models = sorted(models)
+        metafunc.parametrize("model_name", models, ids=models)
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+@pytest.mark.external
+@pytest.mark.parametrize("use_pandas", [True, False], ids=["pandas", "nopandas"])
+def test_load_simulation(function_tmpdir, benchmark, registry, model_name, use_pandas):
+    registry.copy_to(function_tmpdir, model_name)
+    benchmark(lambda: MFSimulation.load(sim_ws=function_tmpdir, use_pandas=use_pandas))
+
+
+@pytest.mark.benchmark
+@pytest.mark.external
+@pytest.mark.slow
+@pytest.mark.parametrize("use_pandas", [True, False], ids=["pandas", "nopandas"])
+def test_write_simulation(function_tmpdir, benchmark, registry, model_name, use_pandas):
+    registry.copy_to(function_tmpdir, model_name)
+    sim = MFSimulation.load(sim_ws=function_tmpdir, use_pandas=use_pandas)
+    benchmark(sim.write_simulation)
diff --git a/autotest/benchmarks/benchmark_mf6listbudget.py b/autotest/benchmarks/benchmark_mf6listbudget.py
new file mode 100644
index 0000000000..7b75ac6845
--- /dev/null
+++ b/autotest/benchmarks/benchmark_mf6listbudget.py
@@ -0,0 +1,33 @@
+import pytest
+
+from flopy.utils.mflistfile import Mf6ListBudget
+
+
+@pytest.fixture
+def mf6_lbf(example_data_path) -> Mf6ListBudget:
+    return Mf6ListBudget(
+        example_data_path / "mf6" / "test001a_Tharmonic" / "flow15.lst"
+    )
+
+
+@pytest.mark.benchmark
+def test_mf6listbudget_load(benchmark, mf6_lbf):
+    benchmark(lambda: Mf6ListBudget(mf6_lbf.file_name))
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize(
+    "incremental", [True, False], ids=["incremental", "no_incremental"]
+)
+def test_mf6listbudget_get_data(benchmark, mf6_lbf, incremental):
+    benchmark(lambda: mf6_lbf.get_data(incremental=incremental))
+
+
+@pytest.mark.benchmark
+def test_mf6listbudget_get_budget(benchmark, mf6_lbf):
+    benchmark(mf6_lbf.get_budget)
+
+
+@pytest.mark.benchmark
+def test_mf6listbudget_to_dataframe(benchmark, mf6_lbf):
+    benchmark(mf6_lbf.get_dataframes)
diff --git a/autotest/benchmarks/benchmark_mflistbudget.py b/autotest/benchmarks/benchmark_mflistbudget.py
new file mode 100644
index 0000000000..4b021c648b
--- /dev/null
+++ b/autotest/benchmarks/benchmark_mflistbudget.py
@@ -0,0 +1,15 @@
+import pytest
+
+from flopy.utils.mflistfile import MfListBudget
+
+
+@pytest.fixture
+def mf_lbf(example_data_path) -> MfListBudget:
+    return MfListBudget(
+        example_data_path / "freyberg_multilayer_transient" / "freyberg.list"
+    )
+
+
+@pytest.mark.benchmark
+def test_mflistbudget_load(benchmark, mf_lbf):
+    benchmark(lambda: MfListBudget(mf_lbf.file_name))
diff --git a/autotest/benchmarks/benchmark_mfusglistbudget.py b/autotest/benchmarks/benchmark_mfusglistbudget.py
new file mode 100644
index 0000000000..0212d62d46
--- /dev/null
+++ b/autotest/benchmarks/benchmark_mfusglistbudget.py
@@ -0,0 +1,19 @@
+import pytest
+
+from flopy.utils.mflistfile import MfusgListBudget
+
+
+@pytest.fixture
+def mfusg_lbf(example_data_path) -> MfusgListBudget:
+    return MfusgListBudget(
+        example_data_path
+        / "mfusg_test"
+        / "03A_conduit_unconfined"
+        / "output"
+        / "ex3A.lst"
+    )
+
+
+@pytest.mark.benchmark
+def test_mfusglistbudget_load(benchmark, mfusg_lbf):
+    benchmark(lambda: MfusgListBudget(mfusg_lbf.file_name))
diff --git a/autotest/benchmarks/benchmark_mtlistfile.py b/autotest/benchmarks/benchmark_mtlistfile.py
new file mode 100644
index 0000000000..a7e87b3f1e
--- /dev/null
+++ b/autotest/benchmarks/benchmark_mtlistfile.py
@@ -0,0 +1,15 @@
+import pytest
+
+from flopy.utils.mtlistfile import MtListBudget
+
+
+@pytest.mark.benchmark
+def test_mtlistfile_load(benchmark, example_data_path):
+    list_file = example_data_path / "mt3d_test" / "mcomp.list"
+
+    def load_and_parse():
+        mt = MtListBudget(str(list_file))
+        mt.parse()
+        return mt.gw_data
+
+    benchmark(load_and_parse)
diff --git a/autotest/benchmarks/benchmark_pathlinefile.py b/autotest/benchmarks/benchmark_pathlinefile.py
new file mode 100644
index 0000000000..d07949f8b3
--- /dev/null
+++ b/autotest/benchmarks/benchmark_pathlinefile.py
@@ -0,0 +1,64 @@
+import pytest
+
+from autotest.test_mp7 import ex01_mf6_model, ex01_mf6_model_name
+from flopy.modpath.mp7 import Modpath7
+from flopy.utils.modpathfile import PathlineFile
+
+
+@pytest.fixture
+def ex01_mp7_model(ex01_mf6_model):
+    sim, function_tmpdir = ex01_mf6_model
+    success, buff = sim.run_simulation()
+    assert success, buff
+    gwf = sim.get_model(ex01_mf6_model_name)
+    mpnam = f"{ex01_mf6_model_name}_mp"
+    mp_ws = function_tmpdir / "mp7"
+    mp_ws.mkdir()
+    return Modpath7.create_mp7(
+        modelname=mpnam,
+        trackdir="forward",
+        flowmodel=gwf,
+        exe_name="mp7",
+        model_ws=mp_ws,
+        rowcelldivisions=1,
+        columncelldivisions=1,
+        layercelldivisions=1,
+    ), mp_ws
+
+
+@pytest.fixture
+def plf(ex01_mp7_model) -> PathlineFile:
+    mp, ws = ex01_mp7_model
+    mp.write_input()
+    success, buff = mp.run_model()
+    assert success
+    return PathlineFile(ws / f"{mp.name}.mppth")
+
+
+@pytest.mark.benchmark
+def test_pathlinefile_load(benchmark, plf):
+    benchmark(lambda: PathlineFile(plf.fname))
+
+
+@pytest.mark.benchmark
+def test_pathlinefile_to_geodataframe(benchmark, ex01_mf6_model, plf):
+    pytest.importorskip("geopandas")
+    sim, function_tmpdir = ex01_mf6_model
+    gwf = sim.get_model()
+    benchmark(lambda: plf.to_geodataframe(gwf.modelgrid))
+
+
+@pytest.mark.benchmark
+def test_pathlinefile_get_data(benchmark, plf):
+    benchmark(plf.get_data)
+
+
+@pytest.mark.benchmark
+def test_pathlinefile_get_alldata(benchmark, plf):
+    benchmark(plf.get_alldata)
+
+
+@pytest.mark.benchmark
+def test_pathlinefile_get_destination_data(benchmark, plf):
+    dest_cells = list(range(0, 100))
+    benchmark(lambda: plf.get_destination_pathline_data(dest_cells))
diff --git a/autotest/benchmarks/benchmark_postprocessing.py b/autotest/benchmarks/benchmark_postprocessing.py
new file mode 100644
index 0000000000..9153c62a63
--- /dev/null
+++ b/autotest/benchmarks/benchmark_postprocessing.py
@@ -0,0 +1,68 @@
+from pathlib import Path
+
+import pytest
+
+from flopy.utils import CellBudgetFile, HeadFile
+from flopy.utils.postprocessing import (
+    get_gradients,
+    get_specific_discharge,
+    get_transmissivities,
+    get_water_table,
+)
+
+from .conftest import load_mf6_sim, load_mf2005_model
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize(
+    "rcxy",
+    [
+        lambda m: (m.dis.nrow.array // 2, m.dis.ncol.array // 2, None, None),
+        lambda m: (
+            None,
+            None,
+            sum(m.modelgrid.extent[:2]) / 2,
+            sum(m.modelgrid.extent[2:4]) / 2,
+        ),
+    ],
+    ids=["r c", "x y"],
+)
+def test_get_transmissivities(benchmark, function_tmpdir, rcxy):
+    sim = load_mf6_sim(function_tmpdir, model_key="freyberg")
+    gwf = sim.get_model()
+    hds_path = Path(function_tmpdir) / "freyberg.hds"
+    hds = HeadFile(hds_path)
+    heads = hds.get_data(totim=hds.get_times()[-1])
+    r, c, x, y = rcxy(gwf)
+    benchmark(lambda: get_transmissivities(heads, gwf, r=r, c=c, x=x, y=y))
+
+
+@pytest.mark.benchmark
+def test_get_water_table(benchmark, function_tmpdir):
+    sim = load_mf6_sim(function_tmpdir, model_key="freyberg")
+    hds_path = Path(function_tmpdir) / "freyberg.hds"
+    hds = HeadFile(hds_path)
+    heads = hds.get_data(totim=hds.get_times()[-1])
+    benchmark(lambda: get_water_table(heads))
+
+
+@pytest.mark.benchmark
+def test_get_gradients(benchmark, function_tmpdir):
+    model = load_mf2005_model(function_tmpdir, model_key="freyberg")
+    hds_path = Path(function_tmpdir) / "freyberg.hds"
+    hds = HeadFile(hds_path)
+    heads = hds.get_data(totim=hds.get_times()[-1])
+    benchmark(lambda: get_gradients(heads, model, nodata=-999))
+
+
+@pytest.mark.benchmark
+def test_get_specific_discharge(benchmark, function_tmpdir):
+    sim = load_mf6_sim(function_tmpdir, model_key="freyberg")
+    gwf = sim.get_model()
+    hds_path = Path(function_tmpdir) / "freyberg.hds"
+    bud_path = Path(function_tmpdir) / "freyberg.cbc"
+    hds = HeadFile(hds_path)
+    bud = CellBudgetFile(bud_path)
+    heads = hds.get_data(totim=hds.get_times()[-1])
+    spdis = bud.get_data(text="DATA-SPDIS")[0]
+    benchmark(lambda: get_specific_discharge(spdis, gwf, heads))
diff --git a/autotest/benchmarks/benchmark_rasters.py b/autotest/benchmarks/benchmark_rasters.py
new file mode 100644
index 0000000000..e8ef8e2d31
--- /dev/null
+++ b/autotest/benchmarks/benchmark_rasters.py
@@ -0,0 +1,136 @@
+import numpy as np
+import pytest
+from modflow_devtools.misc import has_pkg
+
+pytest.importorskip("rasterio")
+pytest.importorskip("affine")
+
+from autotest.conftest import get_example_data_path
+from flopy.discretization.structuredgrid import StructuredGrid
+from flopy.utils.geometry import Polygon
+from flopy.utils.rasters import Raster
+
+
+@pytest.fixture(scope="module")
+def raster_path(example_data_path):
+    return example_data_path / "options" / "dem" / "dem.img"
+
+
+@pytest.fixture(scope="module")
+def raster(raster_path):
+    return Raster.load(raster_path)
+
+
+def origin_and_extent(raster) -> tuple[float, float, float, float]:
+    x0, x1, y0, y1 = raster.bounds
+    # central region
+    cx, cy = (x0 + x1) / 2, (y0 + y1) / 2
+    extent_x, extent_y = (x1 - x0) * 0.6, (y1 - y0) * 0.6
+    grid_x0, grid_y0 = cx - extent_x / 2, cy - extent_y / 2
+    return grid_x0, grid_y0, extent_x, extent_y
+
+
+def make_grid(raster, nrow, ncol) -> StructuredGrid:
+    grid_x0, grid_y0, extent_x, extent_y = origin_and_extent(raster)
+    delr = np.full(ncol, extent_x / ncol)
+    delc = np.full(nrow, extent_y / nrow)
+    return StructuredGrid(
+        delc=delc, delr=delr, xoff=grid_x0, yoff=grid_y0, crs=raster.crs
+    )
+
+
+@pytest.mark.benchmark
+def test_raster_load(benchmark, raster_path):
+    benchmark(lambda: Raster.load(raster_path))
+
+
+RASTER_PATH = get_example_data_path() / "options" / "dem" / "dem.img"
+GRIDS = [
+    make_grid(Raster.load(RASTER_PATH), 10, 10),
+    make_grid(Raster.load(RASTER_PATH), 50, 50),
+    make_grid(Raster.load(RASTER_PATH), 200, 200),
+]
+
+
+@pytest.mark.slow
+@pytest.mark.benchmark
+@pytest.mark.parametrize("grid", GRIDS, ids=["small", "medium", "large"])
+@pytest.mark.parametrize(
+    "method", ["linear", "nearest", "cubic", "mean", "median", "min", "max"]
+)
+def test_raster_resample(benchmark, raster, grid, method):
+    benchmark(lambda: raster.resample_to_grid(grid, band=1, method=method))
+
+
+@pytest.mark.slow
+@pytest.mark.benchmark
+@pytest.mark.skipif(not has_pkg("pyproj"), reason="requires pyproj")
+def test_raster_to_crs_transform(benchmark, raster):
+    benchmark(lambda: raster.to_crs(epsg=4326))
+
+
+def small_poly(raster):
+    # small central polygon, ~20% of extent
+    x0, x1, y0, y1 = raster.bounds
+    cx, cy = (x0 + x1) / 2, (y0 + y1) / 2
+    dx, dy = (x1 - x0) * 0.1, (y1 - y0) * 0.1
+    return Polygon(
+        [(cx - dx, cy - dy), (cx + dx, cy - dy), (cx + dx, cy + dy), (cx - dx, cy + dy)]
+    )
+
+
+def medium_poly(raster):
+    # ~50% of extent
+    x0, x1, y0, y1 = raster.bounds
+    cx, cy = (x0 + x1) / 2, (y0 + y1) / 2
+    dx, dy = (x1 - x0) * 0.25, (y1 - y0) * 0.25
+    return Polygon(
+        [(cx - dx, cy - dy), (cx + dx, cy - dy), (cx + dx, cy + dy), (cx - dx, cy + dy)]
+    )
+
+
+def large_poly(raster):
+    # 80% of extent
+    x0, x1, y0, y1 = raster.bounds
+    dx, dy = (x1 - x0) * 0.1, (y1 - y0) * 0.1
+    return Polygon(
+        [(x0 + dx, y0 + dy), (x1 - dx, y0 + dy), (x1 - dx, y1 - dy), (x0 + dx, y1 - dy)]
+    )
+
+
+POLYGONS = [
+    small_poly(Raster.load(RASTER_PATH)),
+    medium_poly(Raster.load(RASTER_PATH)),
+    large_poly(Raster.load(RASTER_PATH)),
+]
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize("poly", POLYGONS, ids=["small", "medium", "large"])
+def test_raster_crop(benchmark, raster, poly):
+    benchmark(lambda: raster.crop(poly))
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize("poly", POLYGONS, ids=["small", "medium", "large"])
+def test_raster_sample(benchmark, raster, poly):
+    benchmark(lambda: raster.sample_polygon(poly, band=1))
+
+
+@pytest.mark.benchmark
+def test_raster_sample_point(benchmark, raster):
+    x0, x1, y0, y1 = raster.bounds
+    x, y = (x0 + x1) / 2, (y0 + y1) / 2  # center point
+    benchmark(lambda: raster.sample_point(x, y, band=1))
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize("masked", [True, False], ids=["masked", "unmasked"])
+def test_raster_get_array_masked(benchmark, raster, masked):
+    benchmark(lambda: raster.get_array(band=1, masked=masked))
+
+
+@pytest.mark.benchmark
+def test_raster_write(benchmark, raster, function_tmpdir):
+    output_path = function_tmpdir / "output_raster.tif"
+    benchmark(lambda: raster.write(str(output_path)))
diff --git a/autotest/benchmarks/benchmark_sfroutputfile.py b/autotest/benchmarks/benchmark_sfroutputfile.py
new file mode 100644
index 0000000000..7cb5a8c5dc
--- /dev/null
+++ b/autotest/benchmarks/benchmark_sfroutputfile.py
@@ -0,0 +1,35 @@
+import pytest
+
+from flopy.utils.sfroutputfile import SfrFile
+
+
+@pytest.mark.benchmark
+def test_sfrfile_load(benchmark, example_data_path):
+    sfr_file = example_data_path / "sfr_examples" / "test1tr.flw"
+    benchmark(lambda: SfrFile(str(sfr_file)))
+
+
+@pytest.fixture
+def sfrf(example_data_path) -> SfrFile:
+    return SfrFile(str(example_data_path / "sfr_examples" / "test1tr.flw"))
+
+
+@pytest.mark.benchmark
+def test_sfrfile_get_nstrm(benchmark, sfrf):
+    df = sfrf.get_dataframe()
+    benchmark(lambda: SfrFile.get_nstrm(df))
+
+
+@pytest.mark.benchmark
+def test_sfrfile_get_results(benchmark, sfrf):
+    benchmark(lambda: sfrf.get_results(segment=1, reach=1))
+
+
+@pytest.mark.benchmark
+def test_sfrfile_get_times(benchmark, sfrf):
+    benchmark(sfrf.get_times)
+
+
+@pytest.mark.benchmark
+def test_sfrfile_get_dataframe(benchmark, sfrf):
+    benchmark(sfrf.get_dataframe)
diff --git a/autotest/benchmarks/benchmark_ucnfile.py b/autotest/benchmarks/benchmark_ucnfile.py
new file mode 100644
index 0000000000..f1efa814f7
--- /dev/null
+++ b/autotest/benchmarks/benchmark_ucnfile.py
@@ -0,0 +1,32 @@
+import pytest
+
+from flopy.utils import HeadUFile
+
+
+@pytest.mark.benchmark
+def test_headufile_load(benchmark, example_data_path):
+    hds_file = example_data_path / "unstructured" / "headu.githds"
+    benchmark(lambda: HeadUFile(str(hds_file)))
+
+
+@pytest.fixture
+def huf(example_data_path) -> HeadUFile:
+    return HeadUFile(str(example_data_path / "unstructured" / "headu.githds"))
+
+
+@pytest.mark.benchmark
+def test_headufile_get_data(benchmark, huf):
+    times = huf.get_times()
+    mid_time = times[len(times) // 2] if len(times) > 0 else times[0]
+    benchmark(lambda: huf.get_data(totim=mid_time))
+
+
+@pytest.mark.benchmark
+@pytest.mark.slow
+def test_headufile_get_alldata(benchmark, huf):
+    benchmark(huf.get_alldata)
+
+
+@pytest.mark.benchmark
+def test_headufile_get_ts(benchmark, huf):
+    benchmark(lambda: huf.get_ts(0))
diff --git a/autotest/benchmarks/benchmark_zonebudget.py b/autotest/benchmarks/benchmark_zonebudget.py
new file mode 100644
index 0000000000..da0b15b318
--- /dev/null
+++ b/autotest/benchmarks/benchmark_zonebudget.py
@@ -0,0 +1,54 @@
+from pathlib import Path
+
+import numpy as np
+import pytest
+from numpy.typing import NDArray
+
+from flopy.modflow.mf import Modflow
+from flopy.utils.zonbud import ZoneBudget
+
+
+def create_zone_array(nlay, nrow, ncol, n_zones=2) -> NDArray:
+    zones = np.zeros((nlay, nrow, ncol), dtype=np.int32)
+    zone_width = ncol // n_zones  # roughly equal zones
+
+    for i in range(n_zones):
+        start_col = i * zone_width
+        end_col = (i + 1) * zone_width if i < n_zones - 1 else ncol
+        zones[:, :, start_col:end_col] = i + 1
+
+    return zones
+
+
+@pytest.fixture(scope="module", params=[2, 5])
+def case(request, example_data_path) -> tuple[Path, NDArray]:
+    model_path = example_data_path / "zonbud_examples"
+    model = Modflow.load(
+        example_data_path / "freyberg" / "freyberg.nam", version="mf2005"
+    )
+    cbc_path = model_path / "freyberg.gitcbc"
+    zones = create_zone_array(model.nlay, model.nrow, model.ncol, n_zones=request.param)
+    return cbc_path, zones
+
+
+@pytest.mark.slow
+@pytest.mark.benchmark
+def test_zonebudget_load(benchmark, case):
+    cbc_path, zones = case
+    benchmark(lambda: ZoneBudget(str(cbc_path), z=zones))
+
+
+@pytest.mark.slow
+@pytest.mark.benchmark
+def test_zonebudget_get_budget(benchmark, case):
+    cbc_path, zones = case
+    zb = ZoneBudget(str(cbc_path), z=zones)
+    benchmark(zb.get_budget)
+
+
+@pytest.mark.slow
+@pytest.mark.benchmark
+def test_zonebudget_get_dataframes(benchmark, case):
+    cbc_path, zones = case
+    zb = ZoneBudget(str(cbc_path), z=zones)
+    benchmark(zb.get_dataframes)
diff --git a/autotest/benchmarks/conftest.py b/autotest/benchmarks/conftest.py
new file mode 100644
index 0000000000..8ab49b3551
--- /dev/null
+++ b/autotest/benchmarks/conftest.py
@@ -0,0 +1,160 @@
+from pathlib import Path
+from shutil import copytree
+
+import pytest
+from modflow_devtools.misc import run_cmd
+
+from flopy.mf6 import MFSimulation
+from flopy.modflow import Modflow
+
+
+def get_examples_path():
+    """Get path to examples/data directory."""
+    return Path(__file__).parent.parent.parent / "examples" / "data"
+
+
+def load_mf6_sim(tmpdir, model_key="freyberg", run=True):
+    """
+    Load and optionally run a MODFLOW 6 simulation from examples/data.
+
+    Parameters
+    ----------
+    tmpdir : Path
+        Temporary directory to copy the model to
+    model_key : str
+        Model identifier (e.g., "freyberg" for mf6-freyberg)
+    run : bool
+        Whether to run the simulation before returning
+
+    Returns
+    -------
+    MFSimulation
+        The loaded simulation
+    """
+    examples_path = get_examples_path()
+
+    # Map model keys to their actual directory names
+    model_map = {
+        "freyberg": "mf6-freyberg",
+        "test003": "mf6/test003_gwfs_disv",
+        "test006": "mf6/test006_gwf3",
+        "test045": "mf6/test045_lake2tr",
+    }
+
+    model_dir = model_map.get(model_key, f"mf6-{model_key}")
+    source_path = examples_path / model_dir
+
+    if not source_path.exists():
+        raise FileNotFoundError(f"Model directory not found: {source_path}")
+
+    # Copy model files to tmpdir
+    copytree(source_path, tmpdir, dirs_exist_ok=True)
+
+    # Load the simulation
+    sim = MFSimulation.load(sim_ws=tmpdir)
+
+    # Run if requested
+    if run:
+        run_cmd("mf6", cwd=tmpdir)
+
+    return sim
+
+
+def load_mf2005_model(tmpdir, model_key="freyberg", run=False):
+    """
+    Load and optionally run a MODFLOW-2005 model from examples/data.
+
+    Parameters
+    ----------
+    tmpdir : Path
+        Temporary directory to copy the model to
+    model_key : str
+        Model identifier (e.g., "freyberg")
+    run : bool
+        Whether to run the model before returning
+
+    Returns
+    -------
+    Modflow
+        The loaded model
+    """
+    from modflow_devtools.misc import get_namefile_paths
+
+    examples_path = get_examples_path()
+
+    # Map model keys to their actual directory names
+    model_map = {
+        "freyberg": "freyberg_multilayer_transient",
+        "mf2005_test": "mf2005_test",
+    }
+
+    model_dir = model_map.get(model_key, model_key)
+    source_path = examples_path / model_dir
+
+    if not source_path.exists():
+        raise FileNotFoundError(f"Model directory not found: {source_path}")
+
+    # Copy model files to tmpdir
+    copytree(source_path, tmpdir, dirs_exist_ok=True)
+
+    # Find the namefile
+    nam_files = get_namefile_paths(tmpdir, namefile="*.nam")
+    if not nam_files:
+        raise FileNotFoundError(f"No .nam file found in {tmpdir}")
+
+    nam_file = nam_files[0].name
+
+    # Load the model
+    model = Modflow.load(nam_file, model_ws=tmpdir, check=False)
+
+    # Run if requested
+    if run:
+        model.run_model()
+
+    return model
+
+
+@pytest.fixture(scope="session")
+def benchmark_config():
+    """Configure pytest-benchmark settings."""
+    return {
+        "warmup": True,
+        "warmup_iterations": 5,
+        "min_rounds": 5,
+        "disable_gc": True,
+    }
+
+
+@pytest.fixture(scope="session")
+def models_path(request) -> list[Path]:
+    """
+    A directories containing model subdirectories. Use
+    the --models-path command line option once or more to specify
+    model directories. If at least one --models_path is provided,
+    external tests (i.e. those using models from an external repo)
+    will run against model input files found in the given location
+    on the local filesystem rather than model input files from the
+    official model registry. This is useful for testing changes to
+    test model input files during MF6 development.
+    """
+    paths = request.config.getoption("--models-path") or []
+    return [Path(p).expanduser().resolve().absolute() for p in paths]
+
+
+def pytest_addoption(parser):
+    parser.addoption(
+        "--models-path",
+        action="append",
+        type=str,
+        help="directory containing model subdirectories. set this to run external "
+        "tests (i.e. those using models from an external repo) against local model "
+        "input files rather than input files from the official model registry.",
+    )
+    parser.addoption(
+        "--namefile-pattern",
+        action="store",
+        type=str,
+        default="mfsim.nam",
+        help="namefile pattern to use when indexing models when --models-path is set."
+        "does nothing otherwise. default is 'mfsim.nam'.",
+    )
diff --git a/autotest/pytest.ini b/autotest/pytest.ini
deleted file mode 100644
index e867703f14..0000000000
--- a/autotest/pytest.ini
+++ /dev/null
@@ -1,23 +0,0 @@
-[pytest]
-addopts = -ra --color=yes
-python_files =
-    test_*.py
-    profile_*.py
-    benchmark_*.py
-    *_test*.py
-    *_profile*.py
-    *_benchmark*.py
-env_files =
-    .env
-markers =
-    example: exercise scripts, tutorials, notebooks
-    generation: tests for code generation utilities
-    meta: tests run by other tests
-    mf6: tests for MODFLOW 6 support
-    regression: tests comparing multiple versions
-    slow: tests not completing in a few seconds
-filterwarnings =
-    # from python-dateutil, used by arrow, jupyter_client, matplotlib, pandas
-    ignore:datetime.datetime.utcfromtimestamp
-    # from pandas, see https://github.com/pandas-dev/pandas/issues/54466
-    ignore:\n.*Pyarrow
diff --git a/autotest/test_grid_cases.py b/autotest/test_grid_cases.py
index 9415a237bd..24392eedaa 100644
--- a/autotest/test_grid_cases.py
+++ b/autotest/test_grid_cases.py
@@ -58,6 +58,30 @@ def structured_cbd_small():
             laycbd=laycbd,
         )
 
+    @staticmethod
+    def structured_medium():
+        nlay, nrow, ncol = 5, 50, 50
+        return StructuredGrid(
+            delc=np.ones(nrow),
+            delr=np.ones(ncol),
+            top=np.ones((nrow, ncol)) * 100.0,
+            botm=np.array(
+                [np.ones((nrow, ncol)) * (100.0 - (i + 1) * 10.0) for i in range(nlay)]
+            ),
+        )
+
+    @staticmethod
+    def structured_large():
+        nlay, nrow, ncol = 10, 100, 100
+        return StructuredGrid(
+            delc=np.ones(nrow),
+            delr=np.ones(ncol),
+            top=np.ones((nrow, ncol)) * 100.0,
+            botm=np.array(
+                [np.ones((nrow, ncol)) * (100.0 - (i + 1) * 10.0) for i in range(nlay)]
+            ),
+        )
+
     @staticmethod
     def vertex_small():
         nlay, ncpl = 3, 5
diff --git a/pyproject.toml b/pyproject.toml
index 0324b2d054..cce48522c6 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -54,6 +54,7 @@ test = [
     "modflow-devtools>=1.7.0",
     "pytest !=8.1.0",
     "pytest-benchmark",
+    "pytest-codspeed",
     "pytest-cov",
     "pytest-dotenv",
     "pytest-xdist",
@@ -218,3 +219,35 @@ ignore-words-list = [
     "vertx",
     "nd",
 ]
+
+[tool.pytest.ini_options]
+addopts = "-ra --color=yes"
+python_files = [
+    "test_*.py",
+    "profile_*.py",
+    "benchmark_*.py",
+    "*_test*.py",
+    "*_profile*.py",
+    "*_benchmark*.py",
+]
+env_files = [".env"]
+markers = [
+    "benchmark: performance benchmarks",
+    "external: uses models from external repositories",
+    "example: exercise scripts, tutorials, notebooks",
+    "generation: tests for code generation utilities",
+    "meta: tests run by other tests",
+    "mf6: tests for MODFLOW 6 support",
+    "regression: tests comparing multiple versions",
+    "slow: tests not completing in a few seconds",
+]
+filterwarnings = [
+    "ignore:datetime.datetime.utcfromtimestamp",
+    "ignore:\n.*Pyarrow",
+]
+
+[tool.pytest.benchmark]
+# Default settings for pytest-benchmark (ignored by pytest-codspeed)
+min_rounds = 5
+warmup = false
+disable_gc = false