diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index da8ccbf0e..7ed676ce7 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -30,6 +30,7 @@ New features and enhancements
 * Missing values method "pct" and "at_least_n" now accept a new "subfreq" option that allows to compute the missing mask in two steps. When given, the algorithm is applied at this "subfreq" resampling frequency first and then the result is resampled at the target indicator frequency. In the output, a period is invalid if any of its subgroup where flagged as invalid by the chosen method. (:pull:`2058`, :issue:`1820`).
 * ``scipy.stats.rv_continuous`` instances can now be given directly as the ``dist`` argument in ``standardized_precipitation_index`` and ``standardized_precipitation_evapotranspiration_index`` indicators. This includes `lmoments3` distributions when specifying ``method="PWM"``. (:issue:`2043`, :pull:`2045`).
 * Time selection in ``xclim.core.calendar.select_time`` and the ``**indexer`` argument of indicators now supports day-of-year bounds given as DataArrays with spatial and/or temporal dimensions. (:issue:`1987`, :pull:`2055`).
+* Maximum Spacing Estimation method for distribution fitting has been added to `xclim.indices.stats.fit` (:issue:`2078`, :pull:`2077`)
 
 Bug fixes
 ^^^^^^^^^
diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py
index ca2b2874a..9aef2699c 100644
--- a/src/xclim/indices/stats.py
+++ b/src/xclim/indices/stats.py
@@ -55,6 +55,16 @@ def _fitfunc_1d(arr, *, dist, nparams, method, **fitkwargs):
         params = dist.fit(x, *args, method="mle", **kwargs, **fitkwargs)
     elif method == "MM":
         params = dist.fit(x, method="mm", **fitkwargs)
+    elif method in ["MSE", "MPS"]:
+        args, guess = _fit_start(x, dist.name, **fitkwargs)
+        param_info = dist.shapes
+        for i, arg in enumerate(args):
+            guess[param_info[i]] = arg
+
+        fitresult = scipy.stats.fit(
+            dist=dist, data=x, method="mse", guess=guess, **fitkwargs
+        )
+        params = fitresult.params
     elif method == "PWM":
         # lmoments3 will raise an error if only dist.numargs + 2 values are provided
         if len(x) <= dist.numargs + 2:
@@ -107,10 +117,14 @@ def fit(
     dist : str or rv_continuous distribution object
         Name of the univariate distribution, such as beta, expon, genextreme, gamma, gumbel_r, lognorm, norm
         (see :py:mod:scipy.stats for full list) or the distribution object itself.
-    method : {"ML", "MLE", "MM", "PWM", "APP"}
-        Fitting method, either maximum likelihood (ML or MLE), method of moments (MM) or approximate method (APP).
+    method : {"ML", "MLE", "MM", "PWM", "APP", "MSE", "MPS"}
+        Fitting method, either maximum likelihood (ML or MLE), method of moments (MM), maximum product of spacings (MSE or MPS)
+        or approximate method (APP).
         Can also be the probability weighted moments (PWM), also called L-Moments, if a compatible `dist` object is passed.
-        The PWM method is usually more robust to outliers.
+        The PWM method is usually more robust to outliers. The MSE method is more consistent than the MLE method, although
+        it can be more sensitive to repeated data.
+        For the MSE method, each variable parameter must be given finite bounds
+        (provided with keyword argument bounds={'param_name':(min,max),...}).
     dim : str
         The dimension upon which to perform the indexing (default: "time").
     **fitkwargs : dict
@@ -131,6 +145,8 @@ def fit(
         "ML": "maximum likelihood",
         "MM": "method of moments",
         "MLE": "maximum likelihood",
+        "MSE": "maximum product of spacings",
+        "MPS": "maximum product of spacings",
         "PWM": "probability weighted moments",
         "APP": "approximative method",
     }
diff --git a/tests/test_stats.py b/tests/test_stats.py
index 29a40b3f0..96289d717 100644
--- a/tests/test_stats.py
+++ b/tests/test_stats.py
@@ -2,9 +2,12 @@
 
 from __future__ import annotations
 
+from functools import partial
+
 import numpy as np
 import pytest
 import xarray as xr
+from scipy.optimize import differential_evolution
 from scipy.stats import lognorm, norm
 
 from xclim.indices import stats
@@ -157,6 +160,24 @@ def test_genextreme_fit(genextreme):
     np.testing.assert_allclose(p, (0.20949, 297.954091, 75.7911863), 1e-5)
 
 
+def test_mse_fit(genextreme):
+    """Check MSE fit with a series that leads to poor values without good initial conditions."""
+    # Use fixed-seed rng to remove randomness of differential_evolution
+    # (alternative: change optimizer to non-stochastic gradient-based)
+    # TODO: change when minimum scipy exceeds 1.15.0 to rng=0
+    optimizer = partial(differential_evolution, seed=0)
+    p = stats.fit(
+        genextreme,
+        "genextreme",
+        "MSE",
+        bounds=dict(c=(0, 1), scale=(0, 100), loc=(200, 400)),
+        optimizer=optimizer,
+    )
+    np.testing.assert_allclose(
+        p, (0.18435517630019815, 293.61049928703073, 86.70937297745427), 1e-3
+    )
+
+
 def test_fa(fitda):
     T = 10
     q = stats.fa(fitda, T, "lognorm")