diff --git a/doc/user-guide/interpolation.rst b/doc/user-guide/interpolation.rst
index 53d8a52b342..3bc055ae78e 100644
--- a/doc/user-guide/interpolation.rst
+++ b/doc/user-guide/interpolation.rst
@@ -132,10 +132,11 @@ It is now possible to safely compute the difference ``other - interpolated``.
 Interpolation methods
 ---------------------
 
-We use :py:class:`scipy.interpolate.interp1d` for 1-dimensional interpolation.
+We use either :py:class:`scipy.interpolate.interp1d` or special interpolants from
+:py:class:`scipy.interpolate` for 1-dimensional interpolation (see :py:meth:`~xarray.Dataset.interp`).
 For multi-dimensional interpolation, an attempt is first made to decompose the
 interpolation in a series of 1-dimensional interpolations, in which case
-:py:class:`scipy.interpolate.interp1d` is used. If a decomposition cannot be
+the relevant 1-dimensional interpolator is used. If a decomposition cannot be
 made (e.g. with advanced interpolation), :py:func:`scipy.interpolate.interpn` is
 used.
 
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index bb9360b3175..45d27f56f50 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -2228,16 +2228,13 @@ def interp(
         kwargs: Mapping[str, Any] | None = None,
         **coords_kwargs: Any,
     ) -> Self:
-        """Interpolate a DataArray onto new coordinates
+        """
+        Interpolate a DataArray onto new coordinates.
+
+        Performs univariate or multivariate interpolation of a Dataset onto new coordinates,
+        utilizing either NumPy or SciPy interpolation routines.
 
-        Performs univariate or multivariate interpolation of a DataArray onto
-        new coordinates using scipy's interpolation routines. If interpolating
-        along an existing dimension,  either :py:class:`scipy.interpolate.interp1d`
-        or a 1-dimensional scipy interpolator (e.g. :py:class:`scipy.interpolate.KroghInterpolator`)
-        is called. When interpolating along multiple existing dimensions, an
-        attempt is made to decompose the interpolation into multiple
-        1-dimensional interpolations. If this is possible, the 1-dimensional interpolator is called.
-        Otherwise, :py:func:`scipy.interpolate.interpn` is called.
+        Out-of-range values are filled with NaN, unless specified otherwise via `kwargs` to the numpy/scipy interpolant.
 
         Parameters
         ----------
@@ -2246,16 +2243,9 @@ def interp(
             New coordinate can be a scalar, array-like or DataArray.
             If DataArrays are passed as new coordinates, their dimensions are
             used for the broadcasting. Missing values are skipped.
-        method : {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "polynomial"}, default: "linear"
-            The method used to interpolate. The method should be supported by
-            the scipy interpolator:
-
-            - ``interp1d``: {"linear", "nearest", "zero", "slinear",
-              "quadratic", "cubic", "polynomial"}
-            - ``interpn``: {"linear", "nearest"}
-
-            If ``"polynomial"`` is passed, the ``order`` keyword argument must
-            also be provided.
+        method : { "linear", "nearest", "zero", "slinear", "quadratic", "cubic", \
+            "quintic", "polynomial", "pchip", "barycentric", "krogh", "akima", "makima" }
+            Interpolation method to use (see descriptions above).
         assume_sorted : bool, default: False
             If False, values of x can be in any order and they are sorted
             first. If True, x has to be an array of monotonically increasing
@@ -2275,12 +2265,37 @@ def interp(
 
         Notes
         -----
-        scipy is required.
+        - SciPy is required for certain interpolation methods.
+        - When interpolating along multiple dimensions with methods `linear` and `nearest`,
+            the process attempts to decompose the interpolation into independent interpolations
+            along one dimension at a time.
+        - The specific interpolation method and dimensionality determine which
+            interpolant is used:
+
+            1. **Interpolation along one dimension of 1D data (`method='linear'`)**
+                - Uses :py:func:`numpy.interp`, unless `fill_value='extrapolate'` is provided via `kwargs`.
+
+            2. **Interpolation along one dimension of N-dimensional data (N ≥ 1)**
+                - Methods {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "quintic", "polynomial"}
+                    use :py:func:`scipy.interpolate.interp1d`, unless conditions permit the use of :py:func:`numpy.interp`
+                    (as in the case of `method='linear'` for 1D data).
+                - If `method='polynomial'`, the `order` keyword argument must also be provided.
+
+            3. **Special interpolants for interpolation along one dimension of N-dimensional data (N ≥ 1)**
+                - Depending on the `method`, the following interpolants from :py:class:`scipy.interpolate` are used:
+                    - `"pchip"`: :py:class:`scipy.interpolate.PchipInterpolator`
+                    - `"barycentric"`: :py:class:`scipy.interpolate.BarycentricInterpolator`
+                    - `"krogh"`: :py:class:`scipy.interpolate.KroghInterpolator`
+                    - `"akima"` or `"makima"`: :py:class:`scipy.interpolate.Akima1dInterpolator`
+                        (`makima` is handled by passing the `makima` flag).
+
+            4. **Interpolation along multiple dimensions of multi-dimensional data**
+                - Uses :py:func:`scipy.interpolate.interpn` for methods {"linear", "nearest", "slinear",
+                    "cubic", "quintic", "pchip"}.
 
         See Also
         --------
-        scipy.interpolate.interp1d
-        scipy.interpolate.interpn
+        :mod:`scipy.interpolate`
 
         :doc:`xarray-tutorial:fundamentals/02.2_manipulating_dimensions`
             Tutorial material on manipulating data resolution using :py:func:`~xarray.DataArray.interp`
@@ -2376,36 +2391,22 @@ def interp_like(
         """Interpolate this object onto the coordinates of another object,
         filling out of range values with NaN.
 
-        If interpolating along a single existing dimension,
-        :py:class:`scipy.interpolate.interp1d` is called. When interpolating
-        along multiple existing dimensions, an attempt is made to decompose the
-        interpolation into multiple 1-dimensional interpolations. If this is
-        possible, :py:class:`scipy.interpolate.interp1d` is called. Otherwise,
-        :py:func:`scipy.interpolate.interpn` is called.
-
         Parameters
         ----------
         other : Dataset or DataArray
             Object with an 'indexes' attribute giving a mapping from dimension
             names to an 1d array-like, which provides coordinates upon
             which to index the variables in this dataset. Missing values are skipped.
-        method : {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "polynomial"}, default: "linear"
-            The method used to interpolate. The method should be supported by
-            the scipy interpolator:
-
-            - {"linear", "nearest", "zero", "slinear", "quadratic", "cubic",
-              "polynomial"} when ``interp1d`` is called.
-            - {"linear", "nearest"} when ``interpn`` is called.
-
-            If ``"polynomial"`` is passed, the ``order`` keyword argument must
-            also be provided.
+        method : { "linear", "nearest", "zero", "slinear", "quadratic", "cubic", \
+            "quintic", "polynomial", "pchip", "barycentric", "krogh", "akima", "makima" }
+            Interpolation method to use (see descriptions above).
         assume_sorted : bool, default: False
             If False, values of coordinates that are interpolated over can be
             in any order and they are sorted first. If True, interpolated
             coordinates are assumed to be an array of monotonically increasing
             values.
         kwargs : dict, optional
-            Additional keyword passed to scipy's interpolator.
+            Additional keyword arguments passed to the interpolant.
 
         Returns
         -------
@@ -2413,6 +2414,44 @@ def interp_like(
             Another dataarray by interpolating this dataarray's data along the
             coordinates of the other object.
 
+        Notes
+        -----
+        - scipy is required.
+        - If the dataarray has object-type coordinates, reindex is used for these
+            coordinates instead of the interpolation.
+        - When interpolating along multiple dimensions with methods `linear` and `nearest`,
+            the process attempts to decompose the interpolation into independent interpolations
+            along one dimension at a time.
+        - The specific interpolation method and dimensionality determine which
+            interpolant is used:
+
+            1. **Interpolation along one dimension of 1D data (`method='linear'`)**
+                - Uses :py:func:`numpy.interp`, unless `fill_value='extrapolate'` is provided via `kwargs`.
+
+            2. **Interpolation along one dimension of N-dimensional data (N ≥ 1)**
+                - Methods {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "quintic", "polynomial"}
+                    use :py:func:`scipy.interpolate.interp1d`, unless conditions permit the use of :py:func:`numpy.interp`
+                    (as in the case of `method='linear'` for 1D data).
+                - If `method='polynomial'`, the `order` keyword argument must also be provided.
+
+            3. **Special interpolants for interpolation along one dimension of N-dimensional data (N ≥ 1)**
+                - Depending on the `method`, the following interpolants from :py:class:`scipy.interpolate` are used:
+                    - `"pchip"`: :py:class:`scipy.interpolate.PchipInterpolator`
+                    - `"barycentric"`: :py:class:`scipy.interpolate.BarycentricInterpolator`
+                    - `"krogh"`: :py:class:`scipy.interpolate.KroghInterpolator`
+                    - `"akima"` or `"makima"`: :py:class:`scipy.interpolate.Akima1dInterpolator`
+                        (`makima` is handled by passing the `makima` flag).
+
+            4. **Interpolation along multiple dimensions of multi-dimensional data**
+                - Uses :py:func:`scipy.interpolate.interpn` for methods {"linear", "nearest", "slinear",
+                    "cubic", "quintic", "pchip"}.
+
+        See Also
+        --------
+        :func:`DataArray.interp`
+        :func:`DataArray.reindex_like`
+        :mod:`scipy.interpolate`
+
         Examples
         --------
         >>> data = np.arange(12).reshape(4, 3)
@@ -2468,18 +2507,8 @@ def interp_like(
         Coordinates:
           * x        (x) int64 32B 10 20 30 40
           * y        (y) int64 24B 70 80 90
-
-        Notes
-        -----
-        scipy is required.
-        If the dataarray has object-type coordinates, reindex is used for these
-        coordinates instead of the interpolation.
-
-        See Also
-        --------
-        DataArray.interp
-        DataArray.reindex_like
         """
+
         if self.dtype.kind not in "uifc":
             raise TypeError(
                 f"interp only works for a numeric type array. Given {self.dtype}."
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 5341b6f3f6f..3e18c090a2d 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -3909,16 +3909,13 @@ def interp(
         method_non_numeric: str = "nearest",
         **coords_kwargs: Any,
     ) -> Self:
-        """Interpolate a Dataset onto new coordinates
+        """
+        Interpolate a Dataset onto new coordinates.
+
+        Performs univariate or multivariate interpolation of a Dataset onto new coordinates,
+        utilizing either NumPy or SciPy interpolation routines.
 
-        Performs univariate or multivariate interpolation of a Dataset onto
-        new coordinates using scipy's interpolation routines. If interpolating
-        along an existing dimension, either :py:class:`scipy.interpolate.interp1d`
-        or a 1-dimensional scipy interpolator (e.g. :py:class:`scipy.interpolate.KroghInterpolator`)
-        is called.  When interpolating along multiple existing dimensions, an
-        attempt is made to decompose the interpolation into multiple
-        1-dimensional interpolations. If this is possible, the 1-dimensional interpolator
-        is called. Otherwise, :py:func:`scipy.interpolate.interpn` is called.
+        Out-of-range values are filled with NaN, unless specified otherwise via `kwargs` to the numpy/scipy interpolant.
 
         Parameters
         ----------
@@ -3927,28 +3924,17 @@ def interp(
             New coordinate can be a scalar, array-like or DataArray.
             If DataArrays are passed as new coordinates, their dimensions are
             used for the broadcasting. Missing values are skipped.
-        method : {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "polynomial", \
-            "barycentric", "krogh", "pchip", "spline", "akima"}, default: "linear"
-            String indicating which method to use for interpolation:
-
-            - 'linear': linear interpolation. Additional keyword
-              arguments are passed to :py:func:`numpy.interp`
-            - 'nearest', 'zero', 'slinear', 'quadratic', 'cubic', 'polynomial':
-              are passed to :py:func:`scipy.interpolate.interp1d`. If
-              ``method='polynomial'``, the ``order`` keyword argument must also be
-              provided.
-            - 'barycentric', 'krogh', 'pchip', 'spline', 'akima': use their
-              respective :py:class:`scipy.interpolate` classes.
-
+        method : { "linear", "nearest", "zero", "slinear", "quadratic", "cubic", \
+            "quintic", "polynomial", "pchip", "barycentric", "krogh", "akima", "makima" }
+            Interpolation method to use (see descriptions above).
         assume_sorted : bool, default: False
             If False, values of coordinates that are interpolated over can be
             in any order and they are sorted first. If True, interpolated
             coordinates are assumed to be an array of monotonically increasing
             values.
         kwargs : dict, optional
-            Additional keyword arguments passed to scipy's interpolator. Valid
-            options and their behavior depend whether ``interp1d`` or
-            ``interpn`` is used.
+            Additional keyword arguments passed to the interpolator. Valid
+            options and their behavior depend which interpolant is used.
         method_non_numeric : {"nearest", "pad", "ffill", "backfill", "bfill"}, optional
             Method for non-numeric types. Passed on to :py:meth:`Dataset.reindex`.
             ``"nearest"`` is used by default.
@@ -3956,6 +3942,7 @@ def interp(
             The keyword arguments form of ``coords``.
             One of coords or coords_kwargs must be provided.
 
+
         Returns
         -------
         interpolated : Dataset
@@ -3963,12 +3950,37 @@ def interp(
 
         Notes
         -----
-        scipy is required.
+        - SciPy is required for certain interpolation methods.
+        - When interpolating along multiple dimensions with methods `linear` and `nearest`,
+            the process attempts to decompose the interpolation into independent interpolations
+            along one dimension at a time.
+        - The specific interpolation method and dimensionality determine which
+            interpolant is used:
+
+            1. **Interpolation along one dimension of 1D data (`method='linear'`)**
+                - Uses :py:func:`numpy.interp`, unless `fill_value='extrapolate'` is provided via `kwargs`.
+
+            2. **Interpolation along one dimension of N-dimensional data (N ≥ 1)**
+                - Methods {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "quintic", "polynomial"}
+                    use :py:func:`scipy.interpolate.interp1d`, unless conditions permit the use of :py:func:`numpy.interp`
+                    (as in the case of `method='linear'` for 1D data).
+                - If `method='polynomial'`, the `order` keyword argument must also be provided.
+
+            3. **Special interpolants for interpolation along one dimension of N-dimensional data (N ≥ 1)**
+                - Depending on the `method`, the following interpolants from :py:class:`scipy.interpolate` are used:
+                    - `"pchip"`: :py:class:`scipy.interpolate.PchipInterpolator`
+                    - `"barycentric"`: :py:class:`scipy.interpolate.BarycentricInterpolator`
+                    - `"krogh"`: :py:class:`scipy.interpolate.KroghInterpolator`
+                    - `"akima"` or `"makima"`: :py:class:`scipy.interpolate.Akima1dInterpolator`
+                        (`makima` is handled by passing the `makima` flag).
+
+            4. **Interpolation along multiple dimensions of multi-dimensional data**
+                - Uses :py:func:`scipy.interpolate.interpn` for methods {"linear", "nearest", "slinear",
+                    "cubic", "quintic", "pchip"}.
 
         See Also
         --------
-        scipy.interpolate.interp1d
-        scipy.interpolate.interpn
+        :mod:`scipy.interpolate`
 
         :doc:`xarray-tutorial:fundamentals/02.2_manipulating_dimensions`
             Tutorial material on manipulating data resolution using :py:func:`~xarray.Dataset.interp`
@@ -4190,15 +4202,12 @@ def interp_like(
         kwargs: Mapping[str, Any] | None = None,
         method_non_numeric: str = "nearest",
     ) -> Self:
-        """Interpolate this object onto the coordinates of another object,
-        filling the out of range values with NaN.
+        """Interpolate this object onto the coordinates of another object.
 
-        If interpolating along a single existing dimension,
-        :py:class:`scipy.interpolate.interp1d` is called. When interpolating
-        along multiple existing dimensions, an attempt is made to decompose the
-        interpolation into multiple 1-dimensional interpolations. If this is
-        possible, :py:class:`scipy.interpolate.interp1d` is called. Otherwise,
-        :py:func:`scipy.interpolate.interpn` is called.
+        Performs univariate or multivariate interpolation of a Dataset onto new coordinates,
+        utilizing either NumPy or SciPy interpolation routines.
+
+        Out-of-range values are filled with NaN, unless specified otherwise via `kwargs` to the numpy/scipy interpolant.
 
         Parameters
         ----------
@@ -4206,26 +4215,17 @@ def interp_like(
             Object with an 'indexes' attribute giving a mapping from dimension
             names to an 1d array-like, which provides coordinates upon
             which to index the variables in this dataset. Missing values are skipped.
-        method : {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "polynomial", \
-            "barycentric", "krogh", "pchip", "spline", "akima"}, default: "linear"
-            String indicating which method to use for interpolation:
-
-            - 'linear': linear interpolation. Additional keyword
-              arguments are passed to :py:func:`numpy.interp`
-            - 'nearest', 'zero', 'slinear', 'quadratic', 'cubic', 'polynomial':
-              are passed to :py:func:`scipy.interpolate.interp1d`. If
-              ``method='polynomial'``, the ``order`` keyword argument must also be
-              provided.
-            - 'barycentric', 'krogh', 'pchip', 'spline', 'akima': use their
-              respective :py:class:`scipy.interpolate` classes.
-
+        method : { "linear", "nearest", "zero", "slinear", "quadratic", "cubic", \
+            "quintic", "polynomial", "pchip", "barycentric", "krogh", "akima", "makima" }
+            Interpolation method to use (see descriptions above).
         assume_sorted : bool, default: False
             If False, values of coordinates that are interpolated over can be
             in any order and they are sorted first. If True, interpolated
             coordinates are assumed to be an array of monotonically increasing
             values.
         kwargs : dict, optional
-            Additional keyword passed to scipy's interpolator.
+            Additional keyword arguments passed to the interpolator. Valid
+            options and their behavior depend which interpolant is use
         method_non_numeric : {"nearest", "pad", "ffill", "backfill", "bfill"}, optional
             Method for non-numeric types. Passed on to :py:meth:`Dataset.reindex`.
             ``"nearest"`` is used by default.
@@ -4238,14 +4238,41 @@ def interp_like(
 
         Notes
         -----
-        scipy is required.
-        If the dataset has object-type coordinates, reindex is used for these
-        coordinates instead of the interpolation.
+        - scipy is required.
+        - If the dataset has object-type coordinates, reindex is used for these
+            coordinates instead of the interpolation.
+        - When interpolating along multiple dimensions with methods `linear` and `nearest`,
+            the process attempts to decompose the interpolation into independent interpolations
+            along one dimension at a time.
+        - The specific interpolation method and dimensionality determine which
+            interpolant is used:
+
+            1. **Interpolation along one dimension of 1D data (`method='linear'`)**
+                - Uses :py:func:`numpy.interp`, unless `fill_value='extrapolate'` is provided via `kwargs`.
+
+            2. **Interpolation along one dimension of N-dimensional data (N ≥ 1)**
+                - Methods {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "quintic", "polynomial"}
+                    use :py:func:`scipy.interpolate.interp1d`, unless conditions permit the use of :py:func:`numpy.interp`
+                    (as in the case of `method='linear'` for 1D data).
+                - If `method='polynomial'`, the `order` keyword argument must also be provided.
+
+            3. **Special interpolants for interpolation along one dimension of N-dimensional data (N ≥ 1)**
+                - Depending on the `method`, the following interpolants from :py:class:`scipy.interpolate` are used:
+                    - `"pchip"`: :py:class:`scipy.interpolate.PchipInterpolator`
+                    - `"barycentric"`: :py:class:`scipy.interpolate.BarycentricInterpolator`
+                    - `"krogh"`: :py:class:`scipy.interpolate.KroghInterpolator`
+                    - `"akima"` or `"makima"`: :py:class:`scipy.interpolate.Akima1dInterpolator`
+                        (`makima` is handled by passing the `makima` flag).
+
+            4. **Interpolation along multiple dimensions of multi-dimensional data**
+                - Uses :py:func:`scipy.interpolate.interpn` for methods {"linear", "nearest", "slinear",
+                    "cubic", "quintic", "pchip"}.
 
         See Also
         --------
-        Dataset.interp
-        Dataset.reindex_like
+        :func:`Dataset.interp`
+        :func:`Dataset.reindex_like`
+        :mod:`scipy.interpolate`
         """
         if kwargs is None:
             kwargs = {}
diff --git a/xarray/core/missing.py b/xarray/core/missing.py
index 4523e4f8232..ece93cf66ef 100644
--- a/xarray/core/missing.py
+++ b/xarray/core/missing.py
@@ -20,7 +20,7 @@
     timedelta_to_numeric,
 )
 from xarray.core.options import _get_keep_attrs
-from xarray.core.types import Interp1dOptions, InterpOptions
+from xarray.core.types import Interp1dOptions, InterpnOptions, InterpOptions
 from xarray.core.utils import OrderedSet, is_scalar
 from xarray.core.variable import Variable, broadcast_variables
 from xarray.namedarray.parallelcompat import get_chunked_array_type
@@ -154,6 +154,9 @@ def __init__(
                 raise ValueError("order is required when method=polynomial")
             method = order
 
+        if method == "quintic":
+            method = 5
+
         self.method = method
 
         self.cons_kwargs = kwargs
@@ -503,6 +506,8 @@ def _get_interpolator(
             interp_class = _import_interpolant("KroghInterpolator", method)
         elif method == "pchip":
             kwargs.update(axis=-1)
+            # pchip default behavior is to extrapolate
+            kwargs.setdefault("extrapolate", False)
             interp_class = _import_interpolant("PchipInterpolator", method)
         elif method == "spline":
             utils.emit_user_level_warning(
@@ -520,9 +525,6 @@ def _get_interpolator(
         elif method == "makima":
             kwargs.update(method="makima", axis=-1)
             interp_class = _import_interpolant("Akima1DInterpolator", method)
-        elif method == "makima":
-            kwargs.update(method="makima", axis=-1)
-            interp_class = _import_interpolant("Akima1DInterpolator", method)
         else:
             raise ValueError(f"{method} is not a valid scipy interpolator")
     else:
@@ -536,8 +538,7 @@ def _get_interpolator_nd(method, **kwargs):
 
     returns interpolator class and keyword arguments for the class
     """
-    valid_methods = ["linear", "nearest"]
-
+    valid_methods = tuple(get_args(InterpnOptions))
     if method in valid_methods:
         kwargs.update(method=method)
         kwargs.setdefault("bounds_error", False)
@@ -631,8 +632,14 @@ def interp(var, indexes_coords, method: InterpOptions, **kwargs):
         return var.copy()
 
     result = var
-    # decompose the interpolation into a succession of independent interpolation
-    for indep_indexes_coords in decompose_interp(indexes_coords):
+
+    if method in ["linear", "nearest", "slinear"]:
+        # decompose the interpolation into a succession of independent interpolation.
+        indexes_coords = decompose_interp(indexes_coords)
+    else:
+        indexes_coords = [indexes_coords]
+
+    for indep_indexes_coords in indexes_coords:
         var = result
 
         # target dimensions
diff --git a/xarray/core/types.py b/xarray/core/types.py
index 2e7572a3858..3937e4d3631 100644
--- a/xarray/core/types.py
+++ b/xarray/core/types.py
@@ -229,12 +229,20 @@ def copy(
 JoinOptions = Literal["outer", "inner", "left", "right", "exact", "override"]
 
 Interp1dOptions = Literal[
-    "linear", "nearest", "zero", "slinear", "quadratic", "cubic", "polynomial"
+    "linear",
+    "nearest",
+    "zero",
+    "slinear",
+    "quadratic",
+    "cubic",
+    "quintic",
+    "polynomial",
 ]
 InterpolantOptions = Literal[
     "barycentric", "krogh", "pchip", "spline", "akima", "makima"
 ]
-InterpOptions = Union[Interp1dOptions, InterpolantOptions]
+InterpnOptions = Literal["linear", "nearest", "slinear", "cubic", "quintic", "pchip"]
+InterpOptions = Union[Interp1dOptions, InterpolantOptions, InterpnOptions]
 
 DatetimeUnitOptions = Literal[
     "Y", "M", "W", "D", "h", "m", "s", "ms", "us", "μs", "ns", "ps", "fs", "as", None
diff --git a/xarray/tests/test_interp.py b/xarray/tests/test_interp.py
index 6722e8d9404..0e797b73d52 100644
--- a/xarray/tests/test_interp.py
+++ b/xarray/tests/test_interp.py
@@ -1,7 +1,7 @@
 from __future__ import annotations
 
-from itertools import combinations, permutations
-from typing import cast
+from itertools import combinations, permutations, product
+from typing import cast, get_args
 
 import numpy as np
 import pandas as pd
@@ -9,7 +9,12 @@
 
 import xarray as xr
 from xarray.coding.cftimeindex import _parse_array_of_cftime_strings
-from xarray.core.types import InterpOptions
+from xarray.core.types import (
+    Interp1dOptions,
+    InterpnOptions,
+    InterpolantOptions,
+    InterpOptions,
+)
 from xarray.tests import (
     assert_allclose,
     assert_equal,
@@ -28,6 +33,8 @@
 except ImportError:
     pass
 
+ALL_1D = get_args(Interp1dOptions) + get_args(InterpolantOptions)
+
 
 def get_example_data(case: int) -> xr.DataArray:
     if case == 0:
@@ -62,6 +69,43 @@ def get_example_data(case: int) -> xr.DataArray:
         raise ValueError("case must be 1-4")
 
 
+@pytest.fixture
+def nd_interp_coords():
+    # interpolation indices for nd interpolation of da from case 3 of get_example_data
+
+    da = get_example_data(case=3)
+
+    coords = {}
+    # grid -> grid
+    coords["xdestnp"] = np.linspace(0.1, 1.0, 11)
+    coords["ydestnp"] = np.linspace(0.0, 0.2, 10)
+    coords["zdestnp"] = da.z.data
+    # list of the points defined by the above mesh in C order
+    mesh_x, mesh_y, mesh_z = np.meshgrid(
+        coords["xdestnp"], coords["ydestnp"], coords["zdestnp"], indexing="ij"
+    )
+    coords["grid_grid_points"] = np.column_stack(
+        [mesh_x.ravel(), mesh_y.ravel(), mesh_z.ravel()]
+    )
+
+    # grid -> oned
+    coords["xdest"] = xr.DataArray(np.linspace(0.1, 1.0, 11), dims="y")  # type: ignore[assignment]
+    coords["ydest"] = xr.DataArray(np.linspace(0.0, 0.2, 11), dims="y")  # type: ignore[assignment]
+    coords["zdest"] = da.z
+    # grid of the points defined by the oned gridded with zdest in C order
+    coords["grid_oned_points"] = np.array(
+        [
+            (a, b, c)
+            for (a, b), c in product(
+                zip(coords["xdest"].data, coords["ydest"].data, strict=False),
+                coords["zdest"].data,
+            )
+        ]
+    )
+
+    return coords
+
+
 def test_keywargs():
     if not has_scipy:
         pytest.skip("scipy is not installed.")
@@ -245,74 +289,98 @@ def func(obj, dim, new_x, method):
     assert_allclose(actual, expected.transpose("z", "w", "y", transpose_coords=True))
 
 
+@requires_scipy
+@pytest.mark.parametrize("method", get_args(InterpnOptions))
 @pytest.mark.parametrize(
     "case", [pytest.param(3, id="no_chunk"), pytest.param(4, id="chunked")]
 )
-def test_interpolate_nd(case: int) -> None:
-    if not has_scipy:
-        pytest.skip("scipy is not installed.")
-
+def test_interpolate_nd(case: int, method: InterpnOptions, nd_interp_coords) -> None:
     if not has_dask and case == 4:
         pytest.skip("dask is not installed in the environment.")
 
     da = get_example_data(case)
 
     # grid -> grid
-    xdestnp = np.linspace(0.1, 1.0, 11)
-    ydestnp = np.linspace(0.0, 0.2, 10)
-    actual = da.interp(x=xdestnp, y=ydestnp, method="linear")
-
-    # linear interpolation is separateable
-    expected = da.interp(x=xdestnp, method="linear")
-    expected = expected.interp(y=ydestnp, method="linear")
+    xdestnp = nd_interp_coords["xdestnp"]
+    ydestnp = nd_interp_coords["ydestnp"]
+    zdestnp = nd_interp_coords["zdestnp"]
+    grid_grid_points = nd_interp_coords["grid_grid_points"]
+    # the presence/absence of z cordinate may affect nd interpolants, even when the
+    # coordinate is unchanged
+    actual = da.interp(x=xdestnp, y=ydestnp, z=zdestnp, method=method)
+    expected_data = scipy.interpolate.interpn(
+        points=(da.x, da.y, da.z),
+        values=da.load().data,
+        xi=grid_grid_points,
+        method=method,
+        bounds_error=False,
+    ).reshape((len(xdestnp), len(ydestnp), len(zdestnp)))
+    expected = xr.DataArray(
+        expected_data,
+        dims=["x", "y", "z"],
+        coords={
+            "x": xdestnp,
+            "y": ydestnp,
+            "z": zdestnp,
+            "x2": da["x2"].interp(x=xdestnp, method=method),
+        },
+    )
     assert_allclose(actual.transpose("x", "y", "z"), expected.transpose("x", "y", "z"))
 
     # grid -> 1d-sample
-    xdest = xr.DataArray(np.linspace(0.1, 1.0, 11), dims="y")
-    ydest = xr.DataArray(np.linspace(0.0, 0.2, 11), dims="y")
-    actual = da.interp(x=xdest, y=ydest, method="linear")
-
-    # linear interpolation is separateable
-    expected_data = scipy.interpolate.RegularGridInterpolator(
-        (da["x"], da["y"]),
-        da.transpose("x", "y", "z").values,
-        method="linear",
+    xdest = nd_interp_coords["xdest"]
+    ydest = nd_interp_coords["ydest"]
+    zdest = nd_interp_coords["zdest"]
+    grid_oned_points = nd_interp_coords["grid_oned_points"]
+    actual = da.interp(x=xdest, y=ydest, z=zdest, method=method)
+    expected_data = scipy.interpolate.interpn(
+        points=(da.x, da.y, da.z),
+        values=da.data,
+        xi=grid_oned_points,
+        method=method,
         bounds_error=False,
-        fill_value=np.nan,
-    )(np.stack([xdest, ydest], axis=-1))
+    ).reshape([len(xdest), len(zdest)])
     expected = xr.DataArray(
         expected_data,
         dims=["y", "z"],
         coords={
-            "z": da["z"],
             "y": ydest,
+            "z": zdest,
             "x": ("y", xdest.values),
-            "x2": da["x2"].interp(x=xdest),
+            "x2": da["x2"].interp(x=xdest, method=method),
         },
     )
+
     assert_allclose(actual.transpose("y", "z"), expected)
 
     # reversed order
-    actual = da.interp(y=ydest, x=xdest, method="linear")
+    actual = da.interp(y=ydest, x=xdest, z=zdest, method=method)
     assert_allclose(actual.transpose("y", "z"), expected)
 
 
 @requires_scipy
-def test_interpolate_nd_nd() -> None:
+# omit cubic, pchip, quintic because not enough points
+@pytest.mark.parametrize("method", ("linear", "nearest", "slinear"))
+def test_interpolate_nd_nd(method: InterpnOptions) -> None:
     """Interpolate nd array with an nd indexer sharing coordinates."""
     # Create original array
     a = [0, 2]
     x = [0, 1, 2]
-    da = xr.DataArray(
-        np.arange(6).reshape(2, 3), dims=("a", "x"), coords={"a": a, "x": x}
-    )
+    values = np.arange(6).reshape(2, 3)
+    da = xr.DataArray(values, dims=("a", "x"), coords={"a": a, "x": x})
 
     # Create indexer into `a` with dimensions (y, x)
     y = [10]
+    a_targets = [1, 2, 2]
     c = {"x": x, "y": y}
-    ia = xr.DataArray([[1, 2, 2]], dims=("y", "x"), coords=c)
-    out = da.interp(a=ia)
-    expected = xr.DataArray([[1.5, 4, 5]], dims=("y", "x"), coords=c)
+    ia = xr.DataArray([a_targets], dims=("y", "x"), coords=c)
+    out = da.interp(a=ia, method=method)
+
+    expected_xi = list(zip(a_targets, x, strict=False))
+    expected_vec = scipy.interpolate.interpn(
+        points=(a, x), values=values, xi=expected_xi, method=method
+    )
+    expected = xr.DataArray([expected_vec], dims=("y", "x"), coords=c)
     xr.testing.assert_allclose(out.drop_vars("a"), expected)
 
     # If the *shared* indexing coordinates do not match, interp should fail.
@@ -353,14 +421,12 @@ def test_interpolate_nd_with_nan() -> None:
     xr.testing.assert_allclose(out2.drop_vars("a"), expected_ds)
 
 
-@pytest.mark.parametrize("method", ["linear"])
+@requires_scipy
+@pytest.mark.parametrize("method", ("linear",))
 @pytest.mark.parametrize(
     "case", [pytest.param(0, id="no_chunk"), pytest.param(1, id="chunk_y")]
 )
 def test_interpolate_scalar(method: InterpOptions, case: int) -> None:
-    if not has_scipy:
-        pytest.skip("scipy is not installed.")
-
     if not has_dask and case in [1]:
         pytest.skip("dask is not installed in the environment.")
 
@@ -377,6 +443,7 @@ def func(obj, new_x):
             axis=obj.get_axis_num("x"),
             bounds_error=False,
             fill_value=np.nan,
+            kind=method,
         )(new_x)
 
     coords = {"x": xdest, "y": da["y"], "x2": func(da["x2"], xdest)}
@@ -384,33 +451,37 @@ def func(obj, new_x):
     assert_allclose(actual, expected)
 
 
-@pytest.mark.parametrize("method", ["linear"])
+@requires_scipy
+@pytest.mark.parametrize("method", ("linear",))
 @pytest.mark.parametrize(
     "case", [pytest.param(3, id="no_chunk"), pytest.param(4, id="chunked")]
 )
 def test_interpolate_nd_scalar(method: InterpOptions, case: int) -> None:
-    if not has_scipy:
-        pytest.skip("scipy is not installed.")
-
     if not has_dask and case in [4]:
         pytest.skip("dask is not installed in the environment.")
 
     da = get_example_data(case)
     xdest = 0.4
     ydest = 0.05
+    zdest = da.get_index("z")
 
-    actual = da.interp(x=xdest, y=ydest, method=method)
+    actual = da.interp(x=xdest, y=ydest, z=zdest, method=method)
     # scipy interpolation for the reference
     expected_data = scipy.interpolate.RegularGridInterpolator(
-        (da["x"], da["y"]),
+        (da["x"], da["y"], da["z"]),
         da.transpose("x", "y", "z").values,
-        method="linear",
+        method=method,
         bounds_error=False,
         fill_value=np.nan,
-    )(np.stack([xdest, ydest], axis=-1))
-
-    coords = {"x": xdest, "y": ydest, "x2": da["x2"].interp(x=xdest), "z": da["z"]}
-    expected = xr.DataArray(expected_data[0], dims=["z"], coords=coords)
+    )(np.asarray([(xdest, ydest, z_val) for z_val in zdest]))
+
+    coords = {
+        "x": xdest,
+        "y": ydest,
+        "x2": da["x2"].interp(x=xdest, method=method),
+        "z": da["z"],
+    }
+    expected = xr.DataArray(expected_data, dims=["z"], coords=coords)
     assert_allclose(actual, expected)
 
 
@@ -788,7 +859,8 @@ def test_3641() -> None:
 
 
 @requires_scipy
-@pytest.mark.parametrize("method", ["nearest", "linear"])
+# cubic, quintic, pchip omitted because not enough points
+@pytest.mark.parametrize("method", ("linear", "nearest", "slinear"))
 def test_decompose(method: InterpOptions) -> None:
     da = xr.DataArray(
         np.arange(6).reshape(3, 2),
@@ -809,9 +881,7 @@ def test_decompose(method: InterpOptions) -> None:
 
 @requires_scipy
 @requires_dask
-@pytest.mark.parametrize(
-    "method", ["linear", "nearest", "zero", "slinear", "quadratic", "cubic"]
-)
+@pytest.mark.parametrize("method", ("linear", "nearest", "cubic", "pchip", "quintic"))
 @pytest.mark.parametrize("chunked", [True, False])
 @pytest.mark.parametrize(
     "data_ndim,interp_ndim,nscalar",
@@ -830,17 +900,19 @@ def test_interpolate_chunk_1d(
     It should do a series of 1d interpolation
     """
 
+    if method in ["cubic", "pchip", "quintic"] and interp_ndim == 3:
+        pytest.skip("Too slow.")
+
     # 3d non chunked data
-    x = np.linspace(0, 1, 5)
+    x = np.linspace(0, 1, 6)
     y = np.linspace(2, 4, 7)
-    z = np.linspace(-0.5, 0.5, 11)
+    z = np.linspace(-0.5, 0.5, 8)
     da = xr.DataArray(
         data=np.sin(x[:, np.newaxis, np.newaxis])
         * np.cos(y[:, np.newaxis])
         * np.exp(z),
         coords=[("x", x), ("y", y), ("z", z)],
     )
-    kwargs = {"fill_value": "extrapolate"}
 
     # choose the data dimensions
     for data_dims in permutations(da.dims, data_ndim):
@@ -875,8 +947,8 @@ def test_interpolate_chunk_1d(
                         if chunked:
                             dest[dim] = xr.DataArray(data=dest[dim], dims=[dim])
                             dest[dim] = dest[dim].chunk(2)
-                actual = da.interp(method=method, **dest, kwargs=kwargs)
-                expected = da.compute().interp(method=method, **dest, kwargs=kwargs)
+                actual = da.interp(method=method, **dest)
+                expected = da.compute().interp(method=method, **dest)
 
                 assert_identical(actual, expected)
 
@@ -888,7 +960,8 @@ def test_interpolate_chunk_1d(
 
 @requires_scipy
 @requires_dask
-@pytest.mark.parametrize("method", ["linear", "nearest"])
+# quintic omitted because not enough points
+@pytest.mark.parametrize("method", ("linear", "nearest", "slinear", "cubic", "pchip"))
 @pytest.mark.filterwarnings("ignore:Increasing number of chunks")
 def test_interpolate_chunk_advanced(method: InterpOptions) -> None:
     """Interpolate nd array with an nd indexer sharing coordinates."""