diff --git a/CHANGES.rst b/CHANGES.rst
index 1680d3857..bf9636021 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -28,6 +28,7 @@ Bug fixes
 * Fixed and adapted ``time_bnds`` to the newest xarray. (:pull:`1700`).
 * Fixed "agreement fraction" in ``robustness_fractions`` to distinguish between negative change and no change. Added "negative" and "changed negative" fractions (:issue:`1690`, :pull:`1711`).
 * ``make_criteria`` now skips columns with NaNs across all realizations. (:pull:`1713`).
+* Fixed bug QuantileDeltaMapping adjustment not working for seasonal grouping (:issue:`1704`, :pull:`1716`).
 
 Internal changes
 ^^^^^^^^^^^^^^^^
diff --git a/tests/test_sdba/test_adjustment.py b/tests/test_sdba/test_adjustment.py
index 7543ae9d7..3cdc86ff9 100644
--- a/tests/test_sdba/test_adjustment.py
+++ b/tests/test_sdba/test_adjustment.py
@@ -378,6 +378,36 @@ def test_mon_U(
         # Test predict
         np.testing.assert_allclose(p, ref.transpose(*p.dims), rtol=0.1, atol=0.2)
 
+    def test_seasonal(self, series, random):
+        u = random.random(10000)
+        kind = "+"
+        name = "tas"
+        # Define distributions
+        xd = uniform(loc=1, scale=1)
+        yd = uniform(loc=2, scale=4)
+
+        # Generate random numbers with u so we get exact results for comparison
+        x = xd.ppf(u)
+        y = yd.ppf(u)
+
+        # Test train
+        hist = sim = series(x, name)
+        ref = series(y, name)
+
+        QDM = QuantileDeltaMapping.train(
+            ref.astype("float32"),
+            hist.astype("float32"),
+            kind=kind,
+            group="time.season",
+            nquantiles=10,
+        )
+        p = QDM.adjust(sim.astype("float32"), interp="linear")
+
+        # Test predict
+        # Accept discrepancies near extremes
+        middle = (u > 1e-2) * (u < 0.99)
+        np.testing.assert_array_almost_equal(p[middle], ref[middle], 1)
+
     def test_cannon_and_diagnostics(self, cannon_2015_dist, cannon_2015_rvs):
         ref, hist, sim = cannon_2015_rvs(15000, random=False)
 
diff --git a/xclim/sdba/base.py b/xclim/sdba/base.py
index 0d40e716e..6590601ba 100644
--- a/xclim/sdba/base.py
+++ b/xclim/sdba/base.py
@@ -295,6 +295,8 @@ def get_index(
         ind = da.indexes[self.dim]
         if self.prop == "week":
             i = da[self.dim].copy(data=ind.isocalendar().week).astype(int)
+        elif self.prop == "season":
+            i = da[self.dim].copy(data=ind.month % 12 // 3)
         else:
             i = getattr(ind, self.prop)
 
diff --git a/xclim/sdba/utils.py b/xclim/sdba/utils.py
index c36cb247a..09a3129af 100644
--- a/xclim/sdba/utils.py
+++ b/xclim/sdba/utils.py
@@ -360,6 +360,11 @@ def _interp_on_quantiles_2D(newx, newg, oldx, oldy, oldg, method, extrap):  # no
     return out
 
 
+SEASON_MAP = {"DJF": 0, "MAM": 1, "JJA": 2, "SON": 3}
+
+map_season_to_int = np.vectorize(SEASON_MAP.get)
+
+
 @parse_group
 def interp_on_quantiles(
     newx: xr.DataArray,
@@ -429,12 +434,17 @@ def interp_on_quantiles(
         )
         return out
 
-    # else:
     if prop not in xq.dims:
         xq = xq.expand_dims({prop: group.get_coordinate()})
     if prop not in yq.dims:
         yq = yq.expand_dims({prop: group.get_coordinate()})
 
+    # Adding the cyclic bounds fails for string coordinates like seasons
+    # That's why we map the seasons to integers
+    if prop == "season":
+        xq = xq.assign_coords(season=map_season_to_int(xq.season))
+        yq = yq.assign_coords(season=map_season_to_int(yq.season))
+
     xq = add_cyclic_bounds(xq, prop, cyclic_coords=False)
     yq = add_cyclic_bounds(yq, prop, cyclic_coords=False)
     newg = group.get_index(newx, interp=method != "nearest")