diff --git a/NEWS.md b/NEWS.md
index a86a0314..aee2ef07 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,5 @@
+**07/01/2026:** Added `waterdata.parallel_chunks(n)` — a context manager to control how finely the OGC `waterdata` getters split multi-value requests. By default the getters chunk only as much as the server's ~8 KB URL limit forces (the fewest sub-requests). But because every sub-request paginates, splitting a large result further costs little or no extra quota when each sub-request still spans many pages, so a caller who *knows* their pull is large can opt into a finer split for smoother progress, more even concurrency, and a smaller unit of retry/resume: `with waterdata.parallel_chunks(32): df, md = waterdata.get_daily(monitoring_location_id=many_sites)`. `n` is a positive integer (e.g. `2`, `8`, `32`) — the number of sub-requests to fan the call out into; a non-integer or non-positive value raises `ValueError` at the `with`. It caps the *total* sub-request count across every multi-value argument combined (not per argument), bounded below by what the byte limit already forces and above by how many values there are to split, so several multi-value arguments can't multiply past it and `n=1` asks for no extra fan-out. Each sub-request costs a request against your hourly rate limit, and because how many run *at once* is capped separately by `API_USGS_CONCURRENT` (default 32) an `n` beyond that adds quota without adding parallelism, so the useful range is roughly `2` up to `API_USGS_CONCURRENT`. It is deliberately a scoped `with` block rather than an automatic behavior or an environment variable, since the library can't tell in advance whether a query is large — a short-window query might fit in a single page, where extra chunks would only burn quota. Also exported at the top level as `dataretrieval.parallel_chunks`.
+
 **06/23/2026:** **Breaking change (1.2.0):** the minimum supported Python is now **3.10** (`requires-python = ">=3.10"`). 3.9 support was already effectively broken — the `waterdata` module's dependencies (`anyio`, the test stack) require 3.10+, and the `waterdata` test modules already skipped on <3.10. `anyio` is now declared as a direct dependency (it is imported directly by `waterdata`), and the CI/ruff/mypy targets move to 3.10. Also fully removed the deprecated `variable_info` metadata property: the `NWIS_Metadata` override only warned and returned `None` (it relied on the defunct `get_pmcodes`), and the `BaseMetadata` abstract is gone too since nothing implemented it — accessing `.variable_info` now raises `AttributeError`. `site_info` is unaffected.
 
 **06/23/2026:** **Breaking change (1.2.0):** removed the `nadp` module and the deprecated `samples` module ahead of the 1.2.0 release. `nadp` was deprecated on 05/01/2026 — NADP is not a USGS data source, so retrieve NADP data directly from https://nadp.slh.wisc.edu/. The `samples.get_usgs_samples` shim (a deprecated forward to the modern getter) is gone; use `waterdata.get_samples()` instead. `import dataretrieval.nadp` / `import dataretrieval.samples` now raise `ModuleNotFoundError`.
diff --git a/README.md b/README.md
index d904e6d6..73aa273e 100644
--- a/README.md
+++ b/README.md
@@ -105,6 +105,56 @@ df, metadata = waterdata.get_continuous(
 print(f"Retrieved {len(df)} continuous gage height measurements")
 ```
 
+#### Speeding up large downloads with `parallel_chunks`
+
+By default the getters split a multi-value request only as far as the server's
+~8 KB URL limit forces — the fewest sub-requests. For a **large, paginated**
+pull that is needlessly conservative: every sub-request pages through its own
+results, so dividing the query into more, smaller sub-requests lets those pages
+be fetched **in parallel**. `parallel_chunks(n)` opts a single call into that
+finer split, fanning it out into `n` sub-requests. It pays off only when the
+result is large enough to span many pages *and* the query has a multi-value
+argument to divide (such as a list of monitoring locations); on a small query —
+or one with nothing to split — it just adds requests, so it is a deliberate,
+scoped `with` block, never the default.
+
+```python
+from dataretrieval import waterdata
+
+# All stream gages in Ohio, then 20 years of their daily discharge — large
+# enough to span many pages, so it profits from a finer split.
+sites, _ = waterdata.get_monitoring_locations(state="Ohio", site_type_code="ST")
+
+with waterdata.parallel_chunks(32):              # fan out into 32 sub-requests
+    df, md = waterdata.get_daily(
+        monitoring_location_id=sites["monitoring_location_id"].tolist(),
+        parameter_code="00060",                  # discharge
+        time="2004-01-01/2023-12-31",
+    )
+```
+
+`n` is the number of sub-requests to fan the call out into. It is capped by how
+many values there are to split, and each sub-request costs a request against
+your hourly [rate limit](https://api.waterdata.usgs.gov/signup/); since how many
+run *at once* is capped separately by `API_USGS_CONCURRENT` (default 32), the
+useful range is roughly `2` up to that value.
+
+Benchmark — 271 Ohio discharge sites (`get_daily`, `parameter_code="00060"`),
+cold cache, each `n` run against its own time window so results are not
+cache-served, with the page size fixed so every run fetches roughly the same
+number of pages (isolating the effect of parallelism):
+
+| `n`  | parallelism | pages | wall-clock                | speedup |
+| ---- | ----------- | ----- | ------------------------- | ------- |
+| off  | 1           | 32    | 23.9 s / 27.4 s (2 runs)  | 1×      |
+| `8`  | 8           | 37    | 4.1 s / 4.3 s             | ~6×     |
+| `32` | 32          | 44    | 2.0 s                     | ~12×    |
+
+The gain comes from overlapping each sub-request's per-page latency and
+server-side work — a genuinely large pull at the default page size shows a
+similar multi-fold speedup. The extra sub-requests each cost quota, so reserve a
+large `n` for pulls you know are large.
+
 Visit the
 [API Reference](https://doi-usgs.github.io/dataretrieval-python/reference/waterdata.html)
 for more information and examples on available services and input parameters.
diff --git a/dataretrieval/__init__.py b/dataretrieval/__init__.py
index c9df1c45..469fe0f5 100644
--- a/dataretrieval/__init__.py
+++ b/dataretrieval/__init__.py
@@ -44,6 +44,11 @@
     URLTooLong,
 )
 
+# Parallel-chunks control (a context manager). Defined with the chunker in
+# ``dataretrieval.ogc.chunking``; surfaced here for a stable public path
+# ``from dataretrieval import parallel_chunks``.
+from dataretrieval.ogc.chunking import parallel_chunks
+
 # Resumable chunk-interruption exceptions. They are defined in
 # ``dataretrieval.ogc.interruptions`` rather than ``dataretrieval.exceptions``
 # because they carry pandas/httpx state and a resumable ``ChunkedCall`` handle,
@@ -93,5 +98,7 @@
     "ChunkInterrupted",
     "QuotaExhausted",
     "ServiceInterrupted",
+    # parallel-chunks control (defined in ogc.chunking)
+    "parallel_chunks",
     "__version__",
 ]
diff --git a/dataretrieval/ogc/chunking.py b/dataretrieval/ogc/chunking.py
index 00d2766c..c6aa2357 100644
--- a/dataretrieval/ogc/chunking.py
+++ b/dataretrieval/ogc/chunking.py
@@ -9,6 +9,12 @@
 cartesian product of chunks. Requests that already fit get a trivial
 single-step plan — ``ChunkedCall`` has one code path either way.
 
+Parallel chunks: the planner is conservative by default — it splits only as far as
+the byte limit forces. A caller who knows their result is large can opt into a
+finer split via the ``parallel_chunks(n)`` context manager, which fans the query
+out into ``n`` parallel sub-requests; ``n`` drives :meth:`ChunkPlan._refine`. See
+``parallel_chunks`` for the why and the when.
+
 This module owns the *execution* half — the event loop and bounded
 concurrency that drive a plan to completion (``ChunkedCall``) plus the
 public ``multi_value_chunked`` decorator. The neighboring concerns live in
@@ -67,8 +73,10 @@
 
 import asyncio
 import functools
+import numbers
 import os
 from collections.abc import Callable, Iterator
+from contextlib import contextmanager
 from contextvars import copy_context
 from typing import Any, cast
 
@@ -172,6 +180,96 @@ def get_active_client() -> httpx.AsyncClient | None:
     return _chunked_client.get()
 
 
+# Parallel-chunks dial: opt-in to fan a query out *more finely* than the byte
+# limit alone requires. Scoped to a ``with parallel_chunks(...):`` block (a
+# ContextVar), deliberately NOT an env var (see :func:`parallel_chunks` for
+# why). The ambient holds ``n`` — the requested cap on the plan's total
+# sub-request count; ``0`` (the default, outside any block) means "chunk only as
+# much as the byte limit needs".
+_parallel_chunks: Ambient[int] = Ambient("ogc_parallel_chunks", 0)
+
+
+@contextmanager
+def parallel_chunks(n: int) -> Iterator[None]:
+    """
+    Fan the OGC getters' multi-value requests out into ``n`` parallel sub-requests.
+
+    By default the Water Data / NGWMN getters chunk a request only as much as
+    the server's ~8 KB URL-byte limit forces — the fewest sub-requests that
+    fit. That is the safe default, but it can be *needlessly* conservative:
+    because every sub-request paginates, splitting a large result further costs
+    little or no extra quota *as long as each sub-request still spans many
+    pages* — rows-per-chunk far exceeding the page size (ten states pulled as
+    one request then page nearly as many times as ten per-state requests
+    would). When a split leaves each sub-request only a page or two, its partial
+    final page is extra, so finer chunks do add some requests. This context
+    manager lets a caller who *knows* their pull is large ask for that finer
+    split — trading roughly the same pages for more, smaller sub-requests, which
+    gives smoother progress, more even concurrency, and a smaller unit of
+    retry/resume.
+
+    Because the library can't tell in advance whether a query is large (ten
+    states over a short window might fit in a single page, where extra chunks
+    would only burn quota), this is a *deliberate* per-call knob rather than an
+    automatic behavior or a process-wide environment variable — scoping it to a
+    ``with`` block keeps an aggressive setting from leaking into unrelated calls
+    and accidentally spending quota. Outside any block the getters use the
+    conservative default. Only the OGC getters (Water Data, NGWMN) read this;
+    wrapping a legacy NWIS call in the block is a harmless no-op.
+
+    Parameters
+    ----------
+    n : int
+        The number of sub-requests to fan the whole call out into — a positive
+        integer such as ``2``, ``8``, or ``32``. It caps the plan's *total*
+        sub-request count (the cartesian product across every multi-value
+        argument combined, not per argument), so several multi-value arguments
+        cannot multiply past it. The actual count is bounded below by what the
+        ~8 KB URL limit already forces and above by the number of values there
+        are to split, so an ``n`` larger than the input allows simply yields one
+        sub-request per value; ``n=1`` asks for no extra fan-out.
+
+        Each sub-request fetches at least one page, so it costs at least one
+        request against your hourly rate limit — a larger ``n`` spends more
+        quota. And because how many sub-requests run *at once* is capped
+        separately by ``API_USGS_CONCURRENT`` (default 32), an ``n`` beyond that
+        adds quota without adding parallelism; the useful range is roughly ``2``
+        up to ``API_USGS_CONCURRENT``.
+
+    Yields
+    ------
+    None
+
+    Raises
+    ------
+    ValueError
+        If ``n`` is not a positive integer — raised on ``with`` entry, before
+        any request is issued, so a bad value fails loudly rather than silently
+        doing nothing.
+
+    Examples
+    --------
+    >>> from dataretrieval import waterdata
+    >>> with waterdata.parallel_chunks(32):
+    ...     df, md = waterdata.get_daily(
+    ...         monitoring_location_id=many_sites, parameter_code="00060"
+    ...     )  # doctest: +SKIP
+
+    See Also
+    --------
+    ChunkPlan._refine : the planning-side effect of ``n``.
+    """
+    # Accept any integer type (incl. numpy ints, mirroring ``max_rows``); ``bool``
+    # is an ``Integral`` subclass but nonsensical here, so reject it explicitly.
+    if not isinstance(n, numbers.Integral) or isinstance(n, bool) or n < 1:
+        raise ValueError(
+            f"parallel_chunks(n): n must be a positive integer (e.g. 2, 8, 32); "
+            f"got {n!r}."
+        )
+    with _parallel_chunks(n):
+        yield
+
+
 class ChunkedCall:
     """
     Stateful handle for a chunked call.
@@ -591,8 +689,9 @@ def multi_value_chunked(
     ``async def fetch(args) -> (df, response)``, and drives it to
     completion via :meth:`ChunkedCall.resume`. The plan splits multi-value
     list params and the cql-text filter so each sub-request URL fits the
-    byte limit; an already-fitting request is a one-step plan. See the
-    module docstring for the concurrency model.
+    byte limit; an already-fitting request is a one-step plan, unless an
+    active :func:`parallel_chunks` block asks the plan to fan out more
+    finely. See the module docstring for the concurrency model.
 
     Parameters
     ----------
@@ -636,7 +735,14 @@ def wrapper(
             finalize: _Finalize = _passthrough_result,
         ) -> tuple[pd.DataFrame, Any]:
             limit = _OGC_URL_BYTE_LIMIT if url_limit is None else url_limit
-            plan = ChunkPlan(args, build_request, limit)
+            # Read the parallel_chunks dial ``n`` from the ambient set by
+            # ``parallel_chunks`` (0 = off outside any such block; otherwise the
+            # requested total sub-request cap). It only affects *planning*, done
+            # here up front, so a later resume — which re-issues the
+            # already-planned sub-requests — needs no snapshot.
+            plan = ChunkPlan(
+                args, build_request, limit, max_chunks=_parallel_chunks.get()
+            )
             retry_policy = RetryPolicy.from_env()
             # The concurrency cap is resolved inside ``resume()`` from
             # ``API_USGS_CONCURRENT``; ``1`` is a sequential gather,
diff --git a/dataretrieval/ogc/planning.py b/dataretrieval/ogc/planning.py
index 68c337ae..a32e48c6 100644
--- a/dataretrieval/ogc/planning.py
+++ b/dataretrieval/ogc/planning.py
@@ -288,6 +288,21 @@ def _extract_axes(args: dict[str, Any]) -> list[_Axis]:
     return axes
 
 
+def _split_at(chunks: list[list[str]], idx: int) -> None:
+    """Replace ``chunks[idx]`` in place with its two contiguous halves.
+
+    The single primitive both planning passes use to fan an axis out. It
+    preserves the partition invariants every consumer relies on: *coverage*
+    (each atom survives, exactly once) and *contiguous, deterministic order*
+    (resume and :meth:`ChunkPlan.iter_sub_args` depend on it). Kept in one
+    place so those invariants can't drift between :meth:`ChunkPlan._plan`
+    (byte-driven) and :meth:`ChunkPlan._refine` (fan-out-driven).
+    """
+    chunk = chunks[idx]
+    mid = len(chunk) // 2
+    chunks[idx : idx + 1] = [chunk[:mid], chunk[mid:]]
+
+
 class ChunkPlan:
     """
     Strategy for issuing one user-level request as a sequence of
@@ -312,7 +327,18 @@ class ChunkPlan:
         Factory that turns a kwargs dict into a sized httpx request,
         e.g. ``_construct_api_requests``.
     url_limit : int
-        Byte budget for the request (URL + body).
+        Byte budget for the request (URL + body) — a hard ceiling every
+        sub-request must fit.
+    max_chunks : int, optional
+        Soft cap on the plan's total sub-request count (default ``0`` = off).
+        ``0`` chunks only as much as ``url_limit`` requires — the most
+        conservative plan, fewest sub-requests. A positive cap fans the plan
+        out to up to ``max_chunks`` sub-requests overall (the
+        cartesian product across axes, never fewer than the byte budget
+        already forces) — capped as a whole, not per axis, so several
+        multi-value axes can't multiply past the cap. Set from the
+        :func:`~dataretrieval.ogc.chunking.parallel_chunks` ``n``; see
+        :meth:`_refine`.
 
     Attributes
     ----------
@@ -344,6 +370,7 @@ def __init__(
         args: dict[str, Any],
         build_request: Callable[..., httpx.Request],
         url_limit: int,
+        max_chunks: int = 0,
     ) -> None:
         self.args = args
         self.axes: list[_Axis] = []
@@ -352,10 +379,10 @@ def __init__(
 
         axes = _extract_axes(args)
         if not axes:
-            # No chunkable axis: nothing to split. If the single request fits,
-            # run it verbatim (the common passthrough). ``_safe_request_bytes``
-            # treats an un-constructable URL (httpx.InvalidURL, > 64 KB) as over
-            # budget.
+            # No chunkable axis: nothing to split, and ``parallel_chunks`` has
+            # nothing to act on either. If the single request fits, run it
+            # verbatim (the common passthrough). ``_safe_request_bytes`` treats
+            # an un-constructable URL (httpx.InvalidURL, > 64 KB) as over budget.
             if _safe_request_bytes(build_request, args, url_limit) <= url_limit:
                 return
             # Over budget. A filter the chunker doesn't manage — cql-json — is
@@ -388,14 +415,29 @@ def __init__(
         except httpx.InvalidURL:
             initial_request = None
 
+        fits = False
         if initial_request is not None:
             self.canonical_url = str(initial_request.url)
-            if _request_bytes(initial_request) <= url_limit:
-                return
+            fits = _request_bytes(initial_request) <= url_limit
+
+        # A request that already fits and hasn't opted into finer chunking is
+        # the common passthrough: leave ``axes``/``chunks`` empty so
+        # ``total == 1`` and ``iter_sub_args`` yields the original args
+        # verbatim. Only when ``max_chunks`` asks for extra fan-out do
+        # we set the axes up to be refined below.
+        if fits and max_chunks <= 0:
+            return
 
         self.axes = axes
         self.chunks = {axis.arg_key: [list(axis.atoms)] for axis in axes}
-        self._plan(build_request, url_limit)
+        if not fits:
+            # Hard pass: greedy-halve until every worst-case sub-request fits
+            # the byte budget (may raise ``Unchunkable``).
+            self._plan(build_request, url_limit)
+        # Soft pass: optionally split further than the byte budget requires.
+        # Purely additive — never re-raises, and the byte budget stays
+        # satisfied; a no-op at ``max_chunks <= 0``.
+        self._refine(max_chunks)
 
         if self.canonical_url is None:
             # Original URL was un-constructable (httpx.InvalidURL); fall
@@ -447,10 +489,53 @@ def _plan(
                     f"sub-request). Reduce input sizes, shorten or simplify "
                     f"the filter, or split the call manually."
                 )
-            axis_chunks = self.chunks[biggest_axis.arg_key]
-            chunk = axis_chunks[biggest_idx]
-            mid = len(chunk) // 2
-            axis_chunks[biggest_idx : biggest_idx + 1] = [chunk[:mid], chunk[mid:]]
+            _split_at(self.chunks[biggest_axis.arg_key], biggest_idx)
+
+    def _refine(self, max_chunks: int) -> None:
+        """
+        Fan the plan out more finely than the byte budget alone requires —
+        the parallel_chunks dial (see
+        :func:`~dataretrieval.ogc.chunking.parallel_chunks` for why a caller
+        would want this).
+
+        Caps the plan's *total* sub-request count (:attr:`total`, the
+        cartesian product across all axes) at ``max_chunks``, not
+        each axis independently — with several multi-value axes, a cap of 32
+        still means at most 32 sub-requests overall, not ``32 ** n_axes``.
+        Each split picks the single largest splittable chunk across *every*
+        axis (ties broken by axis-extraction order, then lowest index), so
+        growth is distributed round-robin rather than one axis saturating
+        before another is touched. Purely additive — only ever *splits*
+        existing chunks, so the byte pass's work and the ``url_limit``
+        invariant are both preserved, and it never raises. A no-op at
+        ``max_chunks <= 0``.
+
+        Parameters
+        ----------
+        max_chunks : int
+            Soft cap on the plan's total sub-request count (``0`` = off) — the
+            ``parallel_chunks(n)`` value. The cap is on the whole plan (the
+            cartesian product across axes), not per axis, so several multi-value
+            axes can't multiply past it.
+        """
+        if max_chunks <= 0:
+            return
+        while self.total < max_chunks:
+            # Largest splittable chunk across every axis; a chunk of size 1
+            # can't be split further. ``max`` with a stable input order
+            # breaks ties by axis order, then lowest index within an axis.
+            candidate: tuple[_Axis, int] | None = None
+            candidate_size = -1
+            for axis in self.axes:
+                for idx, chunk in enumerate(self.chunks[axis.arg_key]):
+                    if len(chunk) <= 1:
+                        continue
+                    if len(chunk) > candidate_size:
+                        candidate, candidate_size = (axis, idx), len(chunk)
+            if candidate is None:
+                return  # every axis saturated at one atom per chunk
+            axis, idx = candidate
+            _split_at(self.chunks[axis.arg_key], idx)
 
     def _worst_case_args(self) -> dict[str, Any]:
         """
diff --git a/dataretrieval/waterdata/__init__.py b/dataretrieval/waterdata/__init__.py
index 99b6e178..eb231469 100644
--- a/dataretrieval/waterdata/__init__.py
+++ b/dataretrieval/waterdata/__init__.py
@@ -9,6 +9,7 @@
 
 from __future__ import annotations
 
+from dataretrieval.ogc.chunking import parallel_chunks
 from dataretrieval.ogc.filters import FILTER_LANG
 
 # Public API exports
@@ -50,6 +51,7 @@
     "PROFILE_LOOKUP",
     "SERVICES",
     "WATERDATA_SERVICES",
+    "parallel_chunks",
     "get_channel",
     "get_codes",
     "get_combined_metadata",
diff --git a/docs/source/userguide/errors.rst b/docs/source/userguide/errors.rst
index e2dc3ef1..28da515f 100644
--- a/docs/source/userguide/errors.rst
+++ b/docs/source/userguide/errors.rst
@@ -96,6 +96,45 @@ condition clears -- only the unfinished sub-requests are re-issued.
             except ChunkInterrupted as again:
                 exc = again
 
+Chunk a large request more finely
+=================================
+
+By default the getters split an over-large request only as much as the
+server's ~8 KB URL limit forces -- the fewest sub-requests. Because each
+sub-request paginates, splitting a large result further costs little or no
+extra quota *as long as each sub-request still spans many pages* (ten states
+pulled as one request then page nearly as many times as ten per-state requests
+would; a split that leaves each sub-request only a page or two adds its partial
+final page). So if you *know* your pull is large you can ask for a finer split
+with ``parallel_chunks(n)`` -- trading roughly the same pages for more, smaller
+sub-requests, which gives smoother progress, more even concurrency, and a
+smaller unit of retry/resume. It is a scoped ``with``
+block, so an aggressive setting can't leak into unrelated calls and
+accidentally spend quota:
+
+.. code-block:: python
+
+    from dataretrieval import waterdata
+
+    with waterdata.parallel_chunks(32):
+        df, md = waterdata.get_daily(
+            monitoring_location_id=many_sites, parameter_code="00060"
+        )
+
+``n`` is a positive integer (e.g. ``2``, ``8``, ``32``) -- the number of
+sub-requests to fan the call out into; a non-integer or non-positive value
+raises ``ValueError`` at the ``with``. It caps the *total* sub-request count
+across every multi-value argument combined (not per argument), bounded below by
+what the byte limit already forces and above by how many values there are to
+split, so several multi-value arguments can't multiply past it and ``n=1`` asks
+for no extra fan-out. Each sub-request costs a request against your hourly rate
+limit, and because how many run *at once* is capped separately by
+``API_USGS_CONCURRENT`` (default 32) an ``n`` beyond that adds quota without
+adding parallelism -- the useful range is roughly ``2`` up to
+``API_USGS_CONCURRENT``. There is no "off" level: simply don't enter the block
+unless you already expect a large, multi-page result -- on a query that would
+have fit in a single page, extra chunks only burn quota.
+
 The full taxonomy
 =================
 
diff --git a/tests/utils_test.py b/tests/utils_test.py
index b3201419..30950294 100644
--- a/tests/utils_test.py
+++ b/tests/utils_test.py
@@ -208,6 +208,20 @@ def test_chunk_interruptions_exported_at_top_level(self):
             dataretrieval.ChunkInterrupted, dataretrieval.DataRetrievalError
         )
 
+    def test_parallel_chunks_exported_at_top_level_and_waterdata(self):
+        """The ``parallel_chunks`` context manager is reachable both from the top
+        level (``from dataretrieval import parallel_chunks``) and from the
+        user-facing ``dataretrieval.waterdata`` namespace, and both resolve to
+        the single object defined in ``dataretrieval.ogc.chunking``."""
+        import dataretrieval
+        from dataretrieval import waterdata
+        from dataretrieval.ogc import chunking
+
+        assert dataretrieval.parallel_chunks is chunking.parallel_chunks
+        assert waterdata.parallel_chunks is chunking.parallel_chunks
+        assert "parallel_chunks" in dataretrieval.__all__
+        assert "parallel_chunks" in waterdata.__all__
+
 
 class Test_BaseMetadata:
     """Tests of BaseMetadata"""
diff --git a/tests/waterdata_chunking_test.py b/tests/waterdata_chunking_test.py
index 5a946946..850b8709 100644
--- a/tests/waterdata_chunking_test.py
+++ b/tests/waterdata_chunking_test.py
@@ -42,8 +42,10 @@
 from dataretrieval.ogc.chunking import (
     ChunkedCall,
     _chunked_client,
+    _parallel_chunks,
     get_active_client,
     multi_value_chunked,
+    parallel_chunks,
 )
 from dataretrieval.ogc.interruptions import (
     ChunkInterrupted,
@@ -2105,3 +2107,241 @@ async def fetch(args):
     assert "finalized" in df.columns
     assert md[0] == "METADATA"
     assert calls["finalize"] >= 1
+
+
+# ---------------------------------------------------------------------------
+# Parallel chunks: the opt-in dial ``parallel_chunks(n)`` to fan a query out
+# MORE finely than the byte limit alone requires (``ChunkPlan._refine`` + the
+# ``parallel_chunks`` context manager). ``_fake_build``'s base is 200 bytes, so
+# a handful of short atoms sits far under ``url_limit=8000`` — the byte pass
+# passes it through untouched, and any splitting below is the ``n`` cap alone.
+# ``ChunkPlan`` takes the integer cap (``max_chunks``) directly;
+# ``parallel_chunks(n)`` publishes ``n`` onto it. The cap bounds the plan's
+# *total* sub-request count (the cartesian product across axes), not each axis
+# independently — see ``test_cap_caps_the_total_across_axes``.
+# ---------------------------------------------------------------------------
+
+
+def test_zero_cap_preserves_passthrough():
+    """``max_chunks=0`` (the default) must not perturb the existing
+    plan: a multi-value request that fits the byte limit is still the trivial
+    passthrough (no axes, ``total == 1``), byte-for-byte the pre-feature
+    behavior."""
+    args = {"monitoring_location_id": ["A", "B", "C", "D"]}
+    plan = ChunkPlan(args, _fake_build, url_limit=8000, max_chunks=0)
+    assert plan.axes == []
+    assert plan.total == 1
+    assert list(plan.iter_sub_args()) == [args]
+
+
+@pytest.mark.parametrize(
+    ("max_chunks", "expected_pieces"),
+    [(0, 1), (2, 2), (8, 8), (16, 10), (32, 10)],
+)
+def test_cap_ramps_then_saturates(max_chunks, expected_pieces):
+    """A single 10-atom axis that fits the byte limit splits into
+    ``min(10, cap)`` pieces: 1 (off), 2, 8, then saturating at 10 (one atom per
+    chunk) once the cap overshoots the atom count. Monotonic and bounded, and
+    whenever it splits the partition is a cover — every atom exactly once. (The
+    cap-0 passthrough has no axis to cover; see the passthrough test.)"""
+    atoms = [f"S{i:02d}" for i in range(10)]
+    plan = ChunkPlan(
+        {"monitoring_location_id": atoms},
+        _fake_build,
+        url_limit=8000,
+        max_chunks=max_chunks,
+    )
+    assert plan.total == expected_pieces
+    if max_chunks:
+        flattened = [
+            a for chunk in plan.chunks["monitoring_location_id"] for a in chunk
+        ]
+        assert sorted(flattened) == sorted(atoms)
+
+
+def test_cap_bounds_fan_out_for_a_long_axis():
+    """The cap holds fan-out to ``n``: at ``n=32`` a 100-atom axis fans into
+    ``n`` pieces — NOT 100 singletons — so ``parallel_chunks(32)`` on a huge
+    list can't detonate into hundreds of sub-requests. Every atom is still
+    covered exactly once."""
+    high = 32
+    atoms = [f"X{i:03d}" for i in range(100)]
+    plan = ChunkPlan(
+        {"monitoring_location_id": atoms},
+        _fake_build,
+        url_limit=8000,
+        max_chunks=high,
+    )
+    assert plan.total == high
+    flattened = [a for chunk in plan.chunks["monitoring_location_id"] for a in chunk]
+    assert sorted(flattened) == sorted(atoms)
+
+
+def test_cap_below_byte_split_does_not_reduce_fan_out():
+    """The cap is purely additive — it can only split further, never coarsen.
+    A request the byte budget already fans into K>2 chunks is untouched by a
+    cap of 2 (below K), so the byte-driven plan is preserved."""
+    # Heavy axis of four 30-char atoms; a limit tight enough that the byte pass
+    # must drive every atom into its own sub-request (4 pieces > the cap of 2).
+    args = {"monitoring_location_id": ["X" * 30, "Y" * 30, "Z" * 30, "W" * 30]}
+    baseline = ChunkPlan(args, _fake_build, url_limit=250, max_chunks=0)
+    assert baseline.total > 2  # byte pass alone already fanned out past 2
+    refined = ChunkPlan(args, _fake_build, url_limit=250, max_chunks=2)
+    # cap 2 < baseline pieces → refine is a no-op here.
+    assert refined.total == baseline.total
+
+
+def test_cap_never_exceeds_the_byte_budget():
+    """Refining on top of an over-budget request keeps the hard invariant:
+    every sub-request still fits ``url_limit`` (splitting only ever shrinks
+    a chunk), and the fan-out is at least what the byte pass required."""
+    args = {"monitoring_location_id": ["X" * 30, "Y" * 30, "Z" * 30, "W" * 30]}
+    limit = 310
+    byte_only = ChunkPlan(args, _fake_build, url_limit=limit, max_chunks=0)
+    plan = ChunkPlan(args, _fake_build, url_limit=limit, max_chunks=32)
+    assert plan.total >= byte_only.total
+    for sub in plan.iter_sub_args():
+        assert _safe_request_bytes(_fake_build, sub, limit) <= limit
+
+
+def test_cap_refines_the_filter_axis():
+    """The dial treats the cql-text ``filter`` axis like any other: an
+    under-budget filter of N top-level OR-clauses is split along that axis
+    into ``min(N, cap)`` pieces."""
+    clauses = [f"p='{i}'" for i in range(8)]
+    args = {"filter": " OR ".join(clauses)}
+    plan = ChunkPlan(args, _fake_build, url_limit=8000, max_chunks=4)
+    assert len(plan.chunks["filter"]) == 4  # min(8, 4)
+    assert plan.total == 4
+
+
+def test_cap_caps_the_total_across_axes():
+    """With more than one multi-value axis the cap bounds the *total*
+    sub-request count (the cartesian product), not each axis independently —
+    the blast-radius guardrail the dial exists for. Two 6-atom axes at a cap
+    of 4 top out at 4 sub-requests total, not 4x4=16; growth is distributed
+    round-robin across axes rather than one axis alone climbing to the cap."""
+    args = {
+        "monitoring_location_id": [f"L{i}" for i in range(6)],
+        "parameter_code": [f"{i:05d}" for i in range(6)],
+    }
+    plan = ChunkPlan(args, _fake_build, url_limit=8000, max_chunks=4)
+    assert plan.total == 4
+    # Every atom on every axis is still covered exactly once.
+    for key, atoms in (
+        ("monitoring_location_id", args["monitoring_location_id"]),
+        ("parameter_code", args["parameter_code"]),
+    ):
+        flattened = [a for chunk in plan.chunks[key] for a in chunk]
+        assert sorted(flattened) == sorted(atoms)
+
+
+def test_cap_bounds_fan_out_across_many_axes():
+    """The guardrail holds regardless of axis count: three multi-value axes
+    at ``n=32`` still top out at ``n`` sub-requests total, not ``n ** 3`` —
+    the property the single-axis-only cap in the original implementation did
+    not guarantee."""
+    high = 32
+    # Three chunkable axes (two list axes + the filter OR-axis), each with 10
+    # atoms — under the old per-axis cap this would have been high**3.
+    args = {
+        "monitoring_location_id": [f"L{i}" for i in range(10)],
+        "parameter_code": [f"{i:05d}" for i in range(10)],
+        "filter": " OR ".join(f"p='{i}'" for i in range(10)),
+    }
+    plan = ChunkPlan(args, _fake_build, url_limit=8000, max_chunks=high)
+    assert plan.total <= high
+
+
+def test_cap_does_not_mask_unchunkable():
+    """A request with nothing to split that still busts the byte limit must
+    raise ``Unchunkable`` regardless of the cap — the soft pass has no axis to
+    act on and must not swallow the hard failure."""
+    args = {"monitoring_location_id": "one-huge-scalar"}
+    with pytest.raises(Unchunkable):
+        ChunkPlan(args, _fake_build, url_limit=10, max_chunks=32)
+
+
+def test_parallel_chunks_publishes_n_on_the_ambient():
+    """The context manager publishes ``n`` on the ambient for the block and
+    restores the previous value on exit — including proper nesting."""
+    assert _parallel_chunks.get() == 0
+    with parallel_chunks(32):
+        assert _parallel_chunks.get() == 32
+        with parallel_chunks(2):
+            assert _parallel_chunks.get() == 2
+        assert _parallel_chunks.get() == 32  # outer restored
+    assert _parallel_chunks.get() == 0  # default (off) outside any block
+
+
+@pytest.mark.parametrize(
+    "bad",
+    [
+        0,  # not positive
+        -1,  # negative
+        1.5,  # a float, not an int
+        "8",  # a string, even a numeric one
+        "high",  # the old level names are gone
+        None,  # None not accepted
+        True,  # bool is an int subclass but nonsensical here
+        ["8"],  # a list
+    ],
+)
+def test_parallel_chunks_rejects_non_positive_int(bad):
+    """``n`` must be a positive integer; every other shape — zero, negative, a
+    float, a string (including a numeric one and the old level names), ``None``,
+    a ``bool``, a list — raises ``ValueError`` at ``with`` entry, before any
+    request, and leaves the ambient untouched."""
+    with pytest.raises(ValueError, match="must be a positive integer"):
+        with parallel_chunks(bad):
+            pass
+    assert _parallel_chunks.get() == 0
+
+
+def test_parallel_chunks_drives_end_to_end_fan_out():
+    """End-to-end: the same fitting request passes through as a single call by
+    default, but fans into ``n`` sub-requests inside a ``parallel_chunks(n)``
+    block — and the combined result still recovers every atom exactly once."""
+    sites = [f"S{i:02d}" for i in range(8)]
+
+    calls: list[tuple[str, ...]] = []
+
+    @multi_value_chunked(build_request=_fake_build, url_limit=8000)
+    async def fetch(args):
+        chunk = tuple(args["monitoring_location_id"])
+        calls.append(chunk)
+        return pd.DataFrame({"site": list(chunk)}), _ok_response()
+
+    # Default: comfortably under the byte limit → one passthrough call.
+    df_plain, _ = fetch({"monitoring_location_id": sites})
+    assert len(calls) == 1
+    assert sorted(df_plain["site"]) == sorted(sites)
+
+    calls.clear()
+    with parallel_chunks(8):
+        df_fine, _ = fetch({"monitoring_location_id": sites})
+    # 8 atoms at n=8 → 8 singleton sub-requests.
+    assert len(calls) == 8
+    assert all(len(chunk) == 1 for chunk in calls)
+    # Union across chunks recovers the original set, once each.
+    assert sorted(a for chunk in calls for a in chunk) == sorted(sites)
+    assert sorted(df_fine["site"]) == sorted(sites)
+
+
+@pytest.mark.parametrize("n", [1, 2, 3, 8])
+def test_parallel_chunks_supports_arbitrary_n(n):
+    """An arbitrary ``n`` (not only 2/8/32) fans an under-limit request into
+    exactly ``n`` sub-requests, together covering every site once — including
+    ``n=1``, the explicit no-op that stays a single passthrough call."""
+    sites = [f"S{i:02d}" for i in range(8)]
+    calls: list[int] = []
+
+    @multi_value_chunked(build_request=_fake_build, url_limit=8000)
+    async def fetch(args):
+        calls.append(len(args["monitoring_location_id"]))
+        return pd.DataFrame(), _ok_response()
+
+    with parallel_chunks(n):
+        fetch({"monitoring_location_id": sites})
+    assert len(calls) == n
+    assert sum(calls) == 8