diff --git a/CHANGES.rst b/CHANGES.rst
index deebeaeec..7dea66409 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -4,7 +4,7 @@ Changelog
 
 v0.49.0 (unreleased)
 --------------------
-Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user:`juliettelavoie`).
+Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user:`juliettelavoie`), David Huard (:user:`huard`).
 
 Announcements
 ^^^^^^^^^^^^^
@@ -34,6 +34,7 @@ Internal changes
 * Removed the experimental `numba` and `llvm` dependency installation steps in the `tox.ini` file. Added `numba@main` to the upstream dependencies. (:pull:`1709`).
 * Added the `tox-gh` dependency to the development installation recipe. This will soon be required for running the `tox` test ensemble on GitHub Workflows. (:pull:`1709`).
 * Added the `vulture` static code analysis tool for finding dead code to the development dependency list and linters (makefile, tox and pre-commit hooks). (:pull:`1717`).
+* Added error message when using `xclim.indices.stats.dist_method` with `nnlf` and included note in docstring. (:issue:`1683`, :pull:`1714`).
 
 v0.48.2 (2024-02-26)
 --------------------
diff --git a/tests/test_stats.py b/tests/test_stats.py
index adbb6b4c4..18b084fb5 100644
--- a/tests/test_stats.py
+++ b/tests/test_stats.py
@@ -334,7 +334,7 @@ def test_parametric_quantile(use_dask, random):
 
 
 @pytest.mark.parametrize("use_dask", [True, False])
-def test_paramtric_cdf(use_dask, random):
+def test_parametric_cdf(use_dask, random):
     mu = 23
     sigma = 2
     n = 10000
@@ -356,3 +356,16 @@ def test_paramtric_cdf(use_dask, random):
     np.testing.assert_array_almost_equal(out, expected, 1)
     assert "cdf" in out.coords
     assert out.attrs["cell_methods"] == "dparams: cdf"
+
+
+def test_dist_method(fitda):
+    params = stats.fit(fitda, "lognorm")
+    cdf = stats.dist_method(
+        "cdf", fit_params=params, arg=xr.DataArray([0.2, 0.8], dims="val")
+    )
+    assert tuple(cdf.dims) == ("val", "x", "y")
+
+    with pytest.raises(ValueError):
+        stats.dist_method(
+            "nnlf", fit_params=params, dims="val", x=xr.DataArray([0.2, 0.8])
+        )
diff --git a/xclim/indices/stats.py b/xclim/indices/stats.py
index 8c5a559a6..9b72e34e6 100644
--- a/xclim/indices/stats.py
+++ b/xclim/indices/stats.py
@@ -550,8 +550,9 @@ def dist_method(
 ) -> xr.DataArray:
     r"""Vectorized statistical function for given argument on given distribution initialized with params.
 
-    Methods where `"*args"` are the distribution parameters can be wrapped, except those that return new dimensions
-    (Ex: 'rvs' with size != 1, 'stats' with more than one moment, 'interval', 'support')
+    Methods where `"*args"` are the distribution parameters can be wrapped, except those that reduce dimensions (
+    e.g. `nnlf`) or create new dimensions (eg: 'rvs' with size != 1, 'stats' with more than one moment, 'interval',
+    'support').
 
     Parameters
     ----------
@@ -577,7 +578,14 @@ def dist_method(
     """
     # Typically the data to be transformed
     arg = [arg] if arg is not None else []
+    if function == "nnlf":
+        raise ValueError(
+            "This method is not supported because it reduces the dimensionality of the data."
+        )
+
+    # We don't need to set `input_core_dims` because we're explicitly splitting the parameters here.
     args = arg + [fit_params.sel(dparams=dp) for dp in fit_params.dparams.values]
+
     return xr.apply_ufunc(
         _dist_method_1D,
         *args,