diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index b48e5a4041f..3521e8215dd 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -57,6 +57,8 @@ Bug fixes
   `Sam Morley <https://github.com/inakleinbottle>`_.
 - Fixed inconsistencies between docstring and functionality for :py:meth:`DataArray.str.get`
   and :py:meth:`DataArray.str.wrap` (:issue:`4334`). By `Mathias Hauser <https://github.com/mathause>`_.
+- Fixed overflow issue causing incorrect results in computing means of :py:class:`cftime.datetime`
+  arrays (:issue:`4341`). By `Spencer Clark <https://github.com/spencerkclark>`_.
 
 
 Documentation
diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py
index 3d19288228e..377e7377b6a 100644
--- a/xarray/core/duck_array_ops.py
+++ b/xarray/core/duck_array_ops.py
@@ -4,6 +4,7 @@
 accept or return xarray objects.
 """
 import contextlib
+import datetime
 import inspect
 import warnings
 from functools import partial
@@ -470,8 +471,7 @@ def timedelta_to_numeric(value, datetime_unit="ns", dtype=float):
 
 
 def _to_pytimedelta(array, unit="us"):
-    index = pd.TimedeltaIndex(array.ravel(), unit=unit)
-    return index.to_pytimedelta().reshape(array.shape)
+    return array.astype(f"timedelta64[{unit}]").astype(datetime.timedelta)
 
 
 def np_timedelta64_to_float(array, datetime_unit):
diff --git a/xarray/tests/test_duck_array_ops.py b/xarray/tests/test_duck_array_ops.py
index e467629a0b3..7d54aac36f8 100644
--- a/xarray/tests/test_duck_array_ops.py
+++ b/xarray/tests/test_duck_array_ops.py
@@ -333,6 +333,40 @@ def test_cftime_datetime_mean():
     assert_equal(result, expected)
 
 
+@requires_cftime
+def test_cftime_datetime_mean_long_time_period():
+    import cftime
+
+    times = np.array(
+        [
+            [
+                cftime.DatetimeNoLeap(400, 12, 31, 0, 0, 0, 0),
+                cftime.DatetimeNoLeap(520, 12, 31, 0, 0, 0, 0),
+            ],
+            [
+                cftime.DatetimeNoLeap(520, 12, 31, 0, 0, 0, 0),
+                cftime.DatetimeNoLeap(640, 12, 31, 0, 0, 0, 0),
+            ],
+            [
+                cftime.DatetimeNoLeap(640, 12, 31, 0, 0, 0, 0),
+                cftime.DatetimeNoLeap(760, 12, 31, 0, 0, 0, 0),
+            ],
+        ]
+    )
+
+    da = DataArray(times, dims=["time", "d2"])
+    result = da.mean("d2")
+    expected = DataArray(
+        [
+            cftime.DatetimeNoLeap(460, 12, 31, 0, 0, 0, 0),
+            cftime.DatetimeNoLeap(580, 12, 31, 0, 0, 0, 0),
+            cftime.DatetimeNoLeap(700, 12, 31, 0, 0, 0, 0),
+        ],
+        dims=["time"],
+    )
+    assert_equal(result, expected)
+
+
 @requires_cftime
 @requires_dask
 def test_cftime_datetime_mean_dask_error():