diff --git a/xarray/core/_reductions.py b/xarray/core/_reductions.py
index 6b32c60fbaf..240cc655297 100644
--- a/xarray/core/_reductions.py
+++ b/xarray/core/_reductions.py
@@ -17,6 +17,7 @@ def reduce(
         self,
         func: Callable[..., Any],
         dim: Union[None, Hashable, Sequence[Hashable]] = None,
+        *,
         axis: Union[None, int, Sequence[int]] = None,
         keep_attrs: bool = None,
         keepdims: bool = False,
@@ -73,19 +74,8 @@ def count(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.count()
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       int64 5
         """
         return self.reduce(
             duck_array_ops.count,
@@ -144,19 +134,8 @@ def all(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) bool True True True True True False
 
         >>> ds.all()
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       bool False
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -215,19 +194,8 @@ def any(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) bool True True True True True False
 
         >>> ds.any()
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       bool True
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -292,27 +260,12 @@ def max(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.max()
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 3.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.max(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 nan
         """
         return self.reduce(
             duck_array_ops.max,
@@ -378,27 +331,12 @@ def min(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.min()
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 1.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.min(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 nan
         """
         return self.reduce(
             duck_array_ops.min,
@@ -468,27 +406,12 @@ def mean(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.mean()
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 1.8
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.mean(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 nan
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -565,35 +488,16 @@ def prod(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.prod()
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 12.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.prod(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 nan
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.prod(skipna=True, min_count=2)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 12.0
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -671,35 +575,16 @@ def sum(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.sum()
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 9.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.sum(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 nan
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.sum(skipna=True, min_count=2)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 9.0
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -774,35 +659,16 @@ def std(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.std()
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 0.7483
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.std(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 nan
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.std(skipna=True, ddof=1)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 0.8367
         """
         return self.reduce(
             duck_array_ops.std,
@@ -877,35 +743,16 @@ def var(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.var()
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 0.56
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.var(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 nan
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.var(skipna=True, ddof=1)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 0.7
         """
         return self.reduce(
             duck_array_ops.var,
@@ -976,27 +823,12 @@ def median(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.median()
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.median(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  ()
-        Data variables:
-            da       float64 nan
         """
         return self.reduce(
             duck_array_ops.median,
@@ -1015,6 +847,7 @@ def reduce(
         self,
         func: Callable[..., Any],
         dim: Union[None, Hashable, Sequence[Hashable]] = None,
+        *,
         axis: Union[None, int, Sequence[int]] = None,
         keep_attrs: bool = None,
         keepdims: bool = False,
@@ -1070,15 +903,8 @@ def count(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.count()
-        <xarray.DataArray ()>
-        array(5)
         """
         return self.reduce(
             duck_array_ops.count,
@@ -1135,15 +961,8 @@ def all(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ True,  True,  True,  True,  True, False])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.all()
-        <xarray.DataArray ()>
-        array(False)
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -1200,15 +1019,8 @@ def any(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ True,  True,  True,  True,  True, False])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.any()
-        <xarray.DataArray ()>
-        array(True)
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -1271,21 +1083,12 @@ def max(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.max()
-        <xarray.DataArray ()>
-        array(3.)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.max(skipna=False)
-        <xarray.DataArray ()>
-        array(nan)
         """
         return self.reduce(
             duck_array_ops.max,
@@ -1349,21 +1152,12 @@ def min(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.min()
-        <xarray.DataArray ()>
-        array(1.)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.min(skipna=False)
-        <xarray.DataArray ()>
-        array(nan)
         """
         return self.reduce(
             duck_array_ops.min,
@@ -1431,21 +1225,12 @@ def mean(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.mean()
-        <xarray.DataArray ()>
-        array(1.8)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.mean(skipna=False)
-        <xarray.DataArray ()>
-        array(nan)
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -1520,27 +1305,16 @@ def prod(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.prod()
-        <xarray.DataArray ()>
-        array(12.)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.prod(skipna=False)
-        <xarray.DataArray ()>
-        array(nan)
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.prod(skipna=True, min_count=2)
-        <xarray.DataArray ()>
-        array(12.)
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -1616,27 +1390,16 @@ def sum(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.sum()
-        <xarray.DataArray ()>
-        array(9.)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.sum(skipna=False)
-        <xarray.DataArray ()>
-        array(nan)
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.sum(skipna=True, min_count=2)
-        <xarray.DataArray ()>
-        array(9.)
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -1709,27 +1472,16 @@ def std(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.std()
-        <xarray.DataArray ()>
-        array(0.74833148)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.std(skipna=False)
-        <xarray.DataArray ()>
-        array(nan)
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.std(skipna=True, ddof=1)
-        <xarray.DataArray ()>
-        array(0.83666003)
         """
         return self.reduce(
             duck_array_ops.std,
@@ -1802,27 +1554,16 @@ def var(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.var()
-        <xarray.DataArray ()>
-        array(0.56)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.var(skipna=False)
-        <xarray.DataArray ()>
-        array(nan)
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.var(skipna=True, ddof=1)
-        <xarray.DataArray ()>
-        array(0.7)
         """
         return self.reduce(
             duck_array_ops.var,
@@ -1891,21 +1632,12 @@ def median(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.median()
-        <xarray.DataArray ()>
-        array(2.)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.median(skipna=False)
-        <xarray.DataArray ()>
-        array(nan)
         """
         return self.reduce(
             duck_array_ops.median,
@@ -1923,6 +1655,7 @@ def reduce(
         self,
         func: Callable[..., Any],
         dim: Union[None, Hashable, Sequence[Hashable]] = None,
+        *,
         axis: Union[None, int, Sequence[int]] = None,
         keep_attrs: bool = None,
         keepdims: bool = False,
@@ -1979,21 +1712,8 @@ def count(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").count()
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) int64 1 2 2
         """
         return self.reduce(
             duck_array_ops.count,
@@ -2052,21 +1772,8 @@ def all(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) bool True True True True True False
 
         >>> ds.groupby("labels").all()
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) bool False True True
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -2125,21 +1832,8 @@ def any(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) bool True True True True True False
 
         >>> ds.groupby("labels").any()
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) bool True True True
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -2204,31 +1898,12 @@ def max(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").max()
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 1.0 2.0 3.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").max(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 2.0 3.0
         """
         return self.reduce(
             duck_array_ops.max,
@@ -2294,31 +1969,12 @@ def min(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").min()
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 1.0 2.0 1.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").min(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 2.0 1.0
         """
         return self.reduce(
             duck_array_ops.min,
@@ -2388,31 +2044,12 @@ def mean(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").mean()
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 1.0 2.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").mean(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 2.0 2.0
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -2489,41 +2126,16 @@ def prod(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").prod()
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 1.0 4.0 3.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").prod(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 4.0 3.0
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.groupby("labels").prod(skipna=True, min_count=2)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 4.0 3.0
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -2601,41 +2213,16 @@ def sum(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").sum()
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 1.0 4.0 4.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").sum(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 4.0 4.0
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.groupby("labels").sum(skipna=True, min_count=2)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 4.0 4.0
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -2710,41 +2297,16 @@ def std(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").std()
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 0.0 0.0 1.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").std(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 0.0 1.0
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.groupby("labels").std(skipna=True, ddof=1)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 0.0 1.414
         """
         return self.reduce(
             duck_array_ops.std,
@@ -2819,41 +2381,16 @@ def var(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").var()
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 0.0 0.0 1.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").var(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 0.0 1.0
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.groupby("labels").var(skipna=True, ddof=1)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 0.0 2.0
         """
         return self.reduce(
             duck_array_ops.var,
@@ -2924,31 +2461,12 @@ def median(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").median()
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 1.0 2.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").median(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (labels: 3)
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
-        Data variables:
-            da       (labels) float64 nan 2.0 2.0
         """
         return self.reduce(
             duck_array_ops.median,
@@ -2967,6 +2485,7 @@ def reduce(
         self,
         func: Callable[..., Any],
         dim: Union[None, Hashable, Sequence[Hashable]] = None,
+        *,
         axis: Union[None, int, Sequence[int]] = None,
         keep_attrs: bool = None,
         keepdims: bool = False,
@@ -3023,21 +2542,8 @@ def count(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").count()
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) int64 1 3 1
         """
         return self.reduce(
             duck_array_ops.count,
@@ -3096,21 +2602,8 @@ def all(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) bool True True True True True False
 
         >>> ds.resample(time="3M").all()
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) bool True True False
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -3169,21 +2662,8 @@ def any(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) bool True True True True True False
 
         >>> ds.resample(time="3M").any()
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) bool True True True
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -3248,31 +2728,12 @@ def max(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").max()
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 3.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").max(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 3.0 nan
         """
         return self.reduce(
             duck_array_ops.max,
@@ -3338,31 +2799,12 @@ def min(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").min()
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 1.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").min(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 1.0 nan
         """
         return self.reduce(
             duck_array_ops.min,
@@ -3432,31 +2874,12 @@ def mean(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").mean()
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 2.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").mean(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 2.0 nan
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -3533,41 +2956,16 @@ def prod(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").prod()
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 6.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").prod(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 6.0 nan
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.resample(time="3M").prod(skipna=True, min_count=2)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 nan 6.0 nan
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -3645,41 +3043,16 @@ def sum(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").sum()
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 6.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").sum(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 6.0 nan
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.resample(time="3M").sum(skipna=True, min_count=2)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 nan 6.0 nan
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -3754,41 +3127,16 @@ def std(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").std()
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 0.0 0.8165 0.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").std(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 0.0 0.8165 nan
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.resample(time="3M").std(skipna=True, ddof=1)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 nan 1.0 nan
         """
         return self.reduce(
             duck_array_ops.std,
@@ -3863,41 +3211,16 @@ def var(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").var()
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 0.0 0.6667 0.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").var(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 0.0 0.6667 nan
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.resample(time="3M").var(skipna=True, ddof=1)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 nan 1.0 nan
         """
         return self.reduce(
             duck_array_ops.var,
@@ -3968,31 +3291,12 @@ def median(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
-        <xarray.Dataset>
-        Dimensions:  (time: 6)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
-        Data variables:
-            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").median()
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 2.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").median(skipna=False)
-        <xarray.Dataset>
-        Dimensions:  (time: 3)
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
-        Data variables:
-            da       (time) float64 1.0 2.0 nan
         """
         return self.reduce(
             duck_array_ops.median,
@@ -4011,6 +3315,7 @@ def reduce(
         self,
         func: Callable[..., Any],
         dim: Union[None, Hashable, Sequence[Hashable]] = None,
+        *,
         axis: Union[None, int, Sequence[int]] = None,
         keep_attrs: bool = None,
         keepdims: bool = False,
@@ -4066,17 +3371,8 @@ def count(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").count()
-        <xarray.DataArray (labels: 3)>
-        array([1, 2, 2])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.count,
@@ -4133,17 +3429,8 @@ def all(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ True,  True,  True,  True,  True, False])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").all()
-        <xarray.DataArray (labels: 3)>
-        array([False,  True,  True])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -4200,17 +3487,8 @@ def any(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ True,  True,  True,  True,  True, False])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").any()
-        <xarray.DataArray (labels: 3)>
-        array([ True,  True,  True])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -4273,25 +3551,12 @@ def max(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").max()
-        <xarray.DataArray (labels: 3)>
-        array([1., 2., 3.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").max(skipna=False)
-        <xarray.DataArray (labels: 3)>
-        array([nan,  2.,  3.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.max,
@@ -4355,25 +3620,12 @@ def min(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").min()
-        <xarray.DataArray (labels: 3)>
-        array([1., 2., 1.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").min(skipna=False)
-        <xarray.DataArray (labels: 3)>
-        array([nan,  2.,  1.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.min,
@@ -4441,25 +3693,12 @@ def mean(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").mean()
-        <xarray.DataArray (labels: 3)>
-        array([1., 2., 2.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").mean(skipna=False)
-        <xarray.DataArray (labels: 3)>
-        array([nan,  2.,  2.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -4534,33 +3773,16 @@ def prod(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").prod()
-        <xarray.DataArray (labels: 3)>
-        array([1., 4., 3.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").prod(skipna=False)
-        <xarray.DataArray (labels: 3)>
-        array([nan,  4.,  3.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.groupby("labels").prod(skipna=True, min_count=2)
-        <xarray.DataArray (labels: 3)>
-        array([nan,  4.,  3.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -4636,33 +3858,16 @@ def sum(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").sum()
-        <xarray.DataArray (labels: 3)>
-        array([1., 4., 4.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").sum(skipna=False)
-        <xarray.DataArray (labels: 3)>
-        array([nan,  4.,  4.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.groupby("labels").sum(skipna=True, min_count=2)
-        <xarray.DataArray (labels: 3)>
-        array([nan,  4.,  4.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -4735,33 +3940,16 @@ def std(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").std()
-        <xarray.DataArray (labels: 3)>
-        array([0., 0., 1.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").std(skipna=False)
-        <xarray.DataArray (labels: 3)>
-        array([nan,  0.,  1.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.groupby("labels").std(skipna=True, ddof=1)
-        <xarray.DataArray (labels: 3)>
-        array([       nan, 0.        , 1.41421356])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.std,
@@ -4834,33 +4022,16 @@ def var(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").var()
-        <xarray.DataArray (labels: 3)>
-        array([0., 0., 1.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").var(skipna=False)
-        <xarray.DataArray (labels: 3)>
-        array([nan,  0.,  1.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.groupby("labels").var(skipna=True, ddof=1)
-        <xarray.DataArray (labels: 3)>
-        array([nan,  0.,  2.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.var,
@@ -4929,25 +4100,12 @@ def median(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").median()
-        <xarray.DataArray (labels: 3)>
-        array([1., 2., 2.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").median(skipna=False)
-        <xarray.DataArray (labels: 3)>
-        array([nan,  2.,  2.])
-        Coordinates:
-          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.median,
@@ -4965,6 +4123,7 @@ def reduce(
         self,
         func: Callable[..., Any],
         dim: Union[None, Hashable, Sequence[Hashable]] = None,
+        *,
         axis: Union[None, int, Sequence[int]] = None,
         keep_attrs: bool = None,
         keepdims: bool = False,
@@ -5020,17 +4179,8 @@ def count(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").count()
-        <xarray.DataArray (time: 3)>
-        array([1, 3, 1])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.count,
@@ -5087,17 +4237,8 @@ def all(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ True,  True,  True,  True,  True, False])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").all()
-        <xarray.DataArray (time: 3)>
-        array([ True,  True, False])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -5154,17 +4295,8 @@ def any(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ True,  True,  True,  True,  True, False])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").any()
-        <xarray.DataArray (time: 3)>
-        array([ True,  True,  True])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -5227,25 +4359,12 @@ def max(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").max()
-        <xarray.DataArray (time: 3)>
-        array([1., 3., 2.])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").max(skipna=False)
-        <xarray.DataArray (time: 3)>
-        array([ 1.,  3., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.max,
@@ -5309,25 +4428,12 @@ def min(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").min()
-        <xarray.DataArray (time: 3)>
-        array([1., 1., 2.])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").min(skipna=False)
-        <xarray.DataArray (time: 3)>
-        array([ 1.,  1., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.min,
@@ -5395,25 +4501,12 @@ def mean(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").mean()
-        <xarray.DataArray (time: 3)>
-        array([1., 2., 2.])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").mean(skipna=False)
-        <xarray.DataArray (time: 3)>
-        array([ 1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -5488,33 +4581,16 @@ def prod(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").prod()
-        <xarray.DataArray (time: 3)>
-        array([1., 6., 2.])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").prod(skipna=False)
-        <xarray.DataArray (time: 3)>
-        array([ 1.,  6., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.resample(time="3M").prod(skipna=True, min_count=2)
-        <xarray.DataArray (time: 3)>
-        array([nan,  6., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -5590,33 +4666,16 @@ def sum(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").sum()
-        <xarray.DataArray (time: 3)>
-        array([1., 6., 2.])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").sum(skipna=False)
-        <xarray.DataArray (time: 3)>
-        array([ 1.,  6., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.resample(time="3M").sum(skipna=True, min_count=2)
-        <xarray.DataArray (time: 3)>
-        array([nan,  6., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -5689,33 +4748,16 @@ def std(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").std()
-        <xarray.DataArray (time: 3)>
-        array([0.        , 0.81649658, 0.        ])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").std(skipna=False)
-        <xarray.DataArray (time: 3)>
-        array([0.        , 0.81649658,        nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.resample(time="3M").std(skipna=True, ddof=1)
-        <xarray.DataArray (time: 3)>
-        array([nan,  1., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.std,
@@ -5788,33 +4830,16 @@ def var(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").var()
-        <xarray.DataArray (time: 3)>
-        array([0.        , 0.66666667, 0.        ])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").var(skipna=False)
-        <xarray.DataArray (time: 3)>
-        array([0.        , 0.66666667,        nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.resample(time="3M").var(skipna=True, ddof=1)
-        <xarray.DataArray (time: 3)>
-        array([nan,  1., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.var,
@@ -5883,25 +4908,12 @@ def median(
         ...     ),
         ... )
         >>> da
-        <xarray.DataArray (time: 6)>
-        array([ 1.,  2.,  3.,  1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
-            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").median()
-        <xarray.DataArray (time: 3)>
-        array([1., 2., 2.])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").median(skipna=False)
-        <xarray.DataArray (time: 3)>
-        array([ 1.,  2., nan])
-        Coordinates:
-          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.median,
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index ea03ac76f37..1b96f22b744 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -2655,6 +2655,7 @@ def reduce(
         self,
         func: Callable[..., Any],
         dim: Union[None, Hashable, Sequence[Hashable]] = None,
+        *,
         axis: Union[None, int, Sequence[int]] = None,
         keep_attrs: bool = None,
         keepdims: bool = False,
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index e3339f2562f..cf52fed6974 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -5000,6 +5000,7 @@ def reduce(
         self,
         func: Callable,
         dim: Union[Hashable, Iterable[Hashable]] = None,
+        *,
         keep_attrs: bool = None,
         keepdims: bool = False,
         numeric_only: bool = False,
@@ -5035,7 +5036,7 @@ def reduce(
             Dataset with this object's DataArrays replaced with new DataArrays
             of summarized data and the indicated dimension(s) removed.
         """
-        if "axis" in kwargs:
+        if kwargs.get("axis", None) is not None:
             raise ValueError(
                 "passing 'axis' to Dataset reduce methods is ambiguous."
                 " Please use 'dim' instead."
diff --git a/xarray/core/groupby.py b/xarray/core/groupby.py
index 185b4ae5bec..bf5cded33ee 100644
--- a/xarray/core/groupby.py
+++ b/xarray/core/groupby.py
@@ -1,5 +1,6 @@
 import datetime
 import warnings
+from typing import Any, Callable, Hashable, Sequence, Union
 
 import numpy as np
 import pandas as pd
@@ -834,7 +835,15 @@ def _combine(self, applied, shortcut=False):
         return combined
 
     def reduce(
-        self, func, dim=None, axis=None, keep_attrs=None, shortcut=True, **kwargs
+        self,
+        func: Callable[..., Any],
+        dim: Union[None, Hashable, Sequence[Hashable]] = None,
+        axis: Union[None, int, Sequence[int]] = None,
+        keep_attrs: bool = None,
+        keepdims: bool = False,
+        *,
+        shortcut: bool = True,
+        **kwargs: Any,
     ):
         """Reduce the items in this group by applying `func` along some
         dimension(s).
@@ -867,11 +876,15 @@ def reduce(
         if dim is None:
             dim = self._group_dim
 
-        if keep_attrs is None:
-            keep_attrs = _get_keep_attrs(default=False)
-
         def reduce_array(ar):
-            return ar.reduce(func, dim, axis, keep_attrs=keep_attrs, **kwargs)
+            return ar.reduce(
+                func=func,
+                dim=dim,
+                axis=axis,
+                keep_attrs=keep_attrs,
+                keepdims=keepdims,
+                **kwargs,
+            )
 
         check_reduce_dims(dim, self.dims)
 
@@ -949,7 +962,16 @@ def _combine(self, applied):
         combined = self._maybe_unstack(combined)
         return combined
 
-    def reduce(self, func, dim=None, keep_attrs=None, **kwargs):
+    def reduce(
+        self,
+        func: Callable[..., Any],
+        dim: Union[None, Hashable, Sequence[Hashable]] = None,
+        *,
+        axis: Union[None, int, Sequence[int]] = None,
+        keep_attrs: bool = None,
+        keepdims: bool = False,
+        **kwargs: Any,
+    ):
         """Reduce the items in this group by applying `func` along some
         dimension(s).
 
@@ -981,11 +1003,15 @@ def reduce(self, func, dim=None, keep_attrs=None, **kwargs):
         if dim is None:
             dim = self._group_dim
 
-        if keep_attrs is None:
-            keep_attrs = _get_keep_attrs(default=False)
-
         def reduce_dataset(ds):
-            return ds.reduce(func, dim, keep_attrs, **kwargs)
+            return ds.reduce(
+                func=func,
+                dim=dim,
+                axis=axis,
+                keep_attrs=keep_attrs,
+                keepdims=keepdims,
+                **kwargs,
+            )
 
         check_reduce_dims(dim, self.dims)
 
diff --git a/xarray/core/resample.py b/xarray/core/resample.py
index ad11550372c..f1776cf3078 100644
--- a/xarray/core/resample.py
+++ b/xarray/core/resample.py
@@ -1,4 +1,5 @@
 import warnings
+from typing import Any, Callable, Hashable, Sequence, Union
 
 from ._reductions import DataArrayResampleReductions, DatasetResampleReductions
 from .groupby import DataArrayGroupByBase, DatasetGroupByBase
@@ -316,7 +317,15 @@ def apply(self, func, args=(), shortcut=None, **kwargs):
         )
         return self.map(func=func, shortcut=shortcut, args=args, **kwargs)
 
-    def reduce(self, func, dim=None, keep_attrs=None, **kwargs):
+    def reduce(
+        self,
+        func: Callable[..., Any],
+        dim: Union[None, Hashable, Sequence[Hashable]] = None,
+        axis: Union[None, int, Sequence[int]] = None,
+        keep_attrs: bool = None,
+        keepdims: bool = False,
+        **kwargs: Any,
+    ):
         """Reduce the items in this group by applying `func` along the
         pre-defined resampling dimension.
 
@@ -341,4 +350,11 @@ def reduce(self, func, dim=None, keep_attrs=None, **kwargs):
             Array with summarized data and the indicated dimension(s)
             removed.
         """
-        return super().reduce(func, dim, keep_attrs, **kwargs)
+        return super().reduce(
+            func=func,
+            dim=dim,
+            axis=axis,
+            keep_attrs=keep_attrs,
+            keepdims=keepdims,
+            **kwargs,
+        )
diff --git a/xarray/util/generate_reductions.py b/xarray/util/generate_reductions.py
index 4981efcc7e7..fbf2d82dd2c 100644
--- a/xarray/util/generate_reductions.py
+++ b/xarray/util/generate_reductions.py
@@ -38,6 +38,7 @@ def reduce(
         self,
         func: Callable[..., Any],
         dim: Union[None, Hashable, Sequence[Hashable]] = None,
+        *,
         axis: Union[None, int, Sequence[int]] = None,
         keep_attrs: bool = None,
         keepdims: bool = False,