diff --git a/docs/sphinx/source/whatsnew/v0.16.0.rst b/docs/sphinx/source/whatsnew/v0.16.0.rst
index 3ea266e6cb..d7ea94ae8d 100644
--- a/docs/sphinx/source/whatsnew/v0.16.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.16.0.rst
@@ -30,6 +30,10 @@ Breaking Changes
 * Removed the deprecated ``server`` keyword argument from
   :py:func:`pvlib.iotools.sodapro.get_cams`. Use ``url`` instead.
   (:issue:`2767`, :pull:`2766`)
+* Changed the output type of :py:func:`pvlib.iam.marion_ruiz_diffuse`
+  and :py:func:`pvlib.iam.schlick_diffuse` from tuple to ``dict``, to be
+  consistent with :py:func:`pvlib.iam.marion_diffuse`. (:issue:`2837`,
+  :pull:`2842`)
 
 Deprecations
 ~~~~~~~~~~~~
@@ -66,6 +70,9 @@ Enhancements
   (:issue:`2828`, :pull:`2832`)
 * Allow variables from multiple datasets to be requested at once in
   :py:func:`~pvlib.iotools.get_merra2`. (:pull:`2839`)
+* Add support for diffuse IAM in the :py:class:`pvlib.pvsystem.Array` and
+  :py:class:`pvlib.pvsystem.PVSystem` classes (see `pvlib.pvsystem.Array.get_iam_diffuse`
+  and `pvlib.pvsystem.PVSystem.get_iam_diffuse`). (:issue:`2812`, :pull:`2845`)
 
 
 Documentation
diff --git a/pvlib/iam.py b/pvlib/iam.py
index 9ba981c5ea..4c7c4adc91 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -316,9 +316,12 @@ def martin_ruiz(aoi, a_r=0.16):
 
 def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
     '''
-    Determine the incidence angle modifiers (iam) for diffuse sky and
+    Determine the incidence angle modifiers (IAM) for sky diffuse and
     ground-reflected irradiance using the Martin and Ruiz incident angle model.
 
+    As described in [1]_, the IAMs result from integrals that assume the
+    incoming sky diffuse and ground-reflected irradiance are isotropic.
+
     Parameters
     ----------
     surface_tilt: float or array-like, default 0
@@ -346,11 +349,12 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
 
     Returns
     -------
-    iam_sky : numeric
-        The incident angle modifier for sky diffuse
+    iam : dict
+        Incident Angle Modifier (see :term:`iam`) values for each type of
+        diffuse irradiance:
 
-    iam_ground : numeric
-        The incident angle modifier for ground-reflected diffuse
+        * 'sky': radiation from the sky dome
+        * 'ground': radiation reflected from the ground
 
     Notes
     -----
@@ -419,7 +423,9 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
         iam_sky = pd.Series(iam_sky, index=out_index, name='iam_sky')
         iam_gnd = pd.Series(iam_gnd, index=out_index, name='iam_ground')
 
-    return iam_sky, iam_gnd
+    iam = {'sky': iam_sky, 'ground': iam_gnd}
+
+    return iam
 
 
 def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
@@ -583,8 +589,8 @@ def sapm(aoi, module, upper=None):
 
 def marion_diffuse(model, surface_tilt, **kwargs):
     """
-    Determine diffuse irradiance incidence angle modifiers using Marion's
-    method of integrating over solid angle.
+    Determine diffuse irradiance incidence angle modifiers (IAM) using
+    Marion's method of integrating over solid angle.
 
     Parameters
     ----------
@@ -603,7 +609,8 @@ def marion_diffuse(model, surface_tilt, **kwargs):
     Returns
     -------
     iam : dict
-        IAM values for each type of diffuse irradiance:
+        Incident Angle Modifier (see :term:`iam`) values for each type of
+        diffuse irradiance:
 
         * 'sky': radiation from the sky dome (zenith <= 90)
         * 'horizon': radiation from the region of the sky near the horizon
@@ -859,7 +866,7 @@ def schlick(aoi):
 
 def schlick_diffuse(surface_tilt):
     r"""
-    Determine the incidence angle modifiers (IAM) for diffuse sky and
+    Determine the incidence angle modifiers (IAM) for sky diffuse and
     ground-reflected irradiance on a tilted surface using the Schlick
     incident angle model.
 
@@ -890,11 +897,12 @@ def schlick_diffuse(surface_tilt):
 
     Returns
     -------
-    iam_sky : numeric
-        The incident angle modifier for sky diffuse.
+    iam : dict
+        Incident Angle Modifier (see :term:`iam`) values for each type of
+        diffuse irradiance:
 
-    iam_ground : numeric
-        The incident angle modifier for ground-reflected diffuse.
+        * 'sky': radiation from the sky dome
+        * 'ground': radiation reflected from the ground
 
     See Also
     --------
@@ -961,7 +969,9 @@ def schlick_diffuse(surface_tilt):
         cuk = pd.Series(cuk, surface_tilt.index)
         cug = pd.Series(cug, surface_tilt.index)
 
-    return cuk, cug
+    iam = {'sky': cuk, 'ground': cug}
+
+    return iam
 
 
 def _get_model(model_name):
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 9d6ddbe665..73f6ddf5b0 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -697,6 +697,11 @@ def isotropic(surface_tilt, dhi, return_components=False):
             * poa_isotropic: The portion of sky diffuse irradiance on a tilted
               plane from the isotropic sky dome. [Wm⁻²]
 
+    diffuse_components : Dict (array input) or DataFrame (Series input)
+        Keys/columns are:
+            * poa_sky_diffuse: Total sky diffuse
+            * poa_isotropic
+
     References
     ----------
     .. [1] Loutzenhiser P.G. et al. "Empirical validation of models to
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 02b16edaf1..5078f11f15 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -307,7 +307,7 @@ def get_aoi(self, solar_zenith, solar_azimuth):
     @_unwrap_single_value
     def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
                        dni_extra=None, airmass=None, albedo=None,
-                       model='haydavies', **kwargs):
+                       model='haydavies', diffuse_components=False, **kwargs):
         """
         Uses :py:func:`pvlib.irradiance.get_total_irradiance` to
         calculate the plane of array irradiance components on the tilted
@@ -334,6 +334,11 @@ def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
             Ground surface albedo. [unitless]
         model : String, default 'haydavies'
             Irradiance model.
+        diffuse_components : bool, default False
+            If `True`, returns values for the different diffuse irradiance
+            components available from the selected model
+            (e.g., isotropic, circumsolar, horizon brightening).
+            If `False`, only the total diffuse irradiance is returned.
 
         kwargs
             Extra parameters passed to
@@ -373,7 +378,9 @@ def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
             array.get_irradiance(solar_zenith, solar_azimuth,
                                  dni, ghi, dhi,
                                  dni_extra=dni_extra, airmass=airmass,
-                                 albedo=albedo, model=model, **kwargs)
+                                 albedo=albedo, model=model,
+                                 diffuse_components=diffuse_components,
+                                 **kwargs)
             for array, dni, ghi, dhi, albedo in zip(
                 self.arrays, dni, ghi, dhi, albedo
             )
@@ -382,8 +389,8 @@ def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
     @_unwrap_single_value
     def get_iam(self, aoi, iam_model='physical'):
         """
-        Determine the incidence angle modifier using the method specified by
-        ``iam_model``.
+        Determine the incidence angle modifier for direct irradiance
+        using the method specified by ``iam_model``.
 
         Parameters for the selected IAM model are expected to be in
         ``PVSystem.module_parameters``. Default parameters are available for
@@ -411,6 +418,48 @@ def get_iam(self, aoi, iam_model='physical'):
         return tuple(array.get_iam(aoi, iam_model)
                      for array, aoi in zip(self.arrays, aoi))
 
+    @_unwrap_single_value
+    def get_iam_diffuse(self, surface_tilt, iam_model='marion_diffuse',
+                        marion_model=None, **kwargs):
+        """
+        Determine the incidence angle modifier for diffuse irradiance using the
+        method specified by ``iam_model``.
+
+        Parameters for the selected IAM model are expected to be in
+        ``Array.module_parameters``. Default parameters are available for
+        the 'marion_diffuse' and 'martin_ruiz_diffuse' models.
+
+        Parameters
+        ----------
+        surface_tilt : float or Series
+            The tilt angle of the surface in degrees.
+        iam_model : string, default 'marion_diffuse'
+            The IAM model to be used. Valid strings are 'marion_diffuse',
+            'martin_ruiz_diffuse', and 'schlick_diffuse'.
+        marion_model : string, default None
+            The IAM function to evaluate across a solid angle. Only used when
+            ``iam_model='marion_diffuse'``. Must be one of 'ashrae',
+            'physical', 'martin_ruiz', 'sapm', and 'schlick'.
+
+        kwargs : dict, optional
+            Additional keyword arguments passed to the IAM model function.
+
+        Returns
+        -------
+        iam_diffuse : dict or DataFrame
+            The AOI modifiers for different diffuse irradiance components.
+            Included components depend on the selected ``iam_model``.
+
+        Raises
+        ------
+        ValueError
+            if `iam_model` is not a valid model name.
+        """
+        surface_tilt = self._validate_per_array(surface_tilt)
+        return tuple(array.get_iam_diffuse(tilt, iam_model=iam_model,
+                                           marion_model=marion_model, **kwargs)
+                     for array, tilt in zip(self.arrays, surface_tilt))
+
     @_unwrap_single_value
     def get_cell_temperature(self, poa_global, temp_air, wind_speed, model,
                              effective_irradiance=None, longwave_down=None):
@@ -1093,7 +1142,7 @@ def get_aoi(self, solar_zenith, solar_azimuth):
 
     def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
                        dni_extra=None, airmass=None, albedo=None,
-                       model='haydavies', **kwargs):
+                       model='haydavies', diffuse_components=False, **kwargs):
         """
         Get plane of array irradiance components.
 
@@ -1121,6 +1170,11 @@ def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
             Ground surface albedo. [unitless]
         model : String, default 'haydavies'
             Irradiance model.
+        diffuse_components : bool, default False
+            If `True`, returns values for the different diffuse irradiance
+            components available from the selected model
+            (e.g., isotropic, circumsolar, horizon brightening).
+            If `False`, only the total diffuse irradiance is returned.
 
         kwargs
             Extra parameters passed to
@@ -1161,20 +1215,23 @@ def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
             airmass = atmosphere.get_relative_airmass(solar_zenith)
 
         orientation = self.mount.get_orientation(solar_zenith, solar_azimuth)
-        return irradiance.get_total_irradiance(orientation['surface_tilt'],
-                                               orientation['surface_azimuth'],
-                                               solar_zenith, solar_azimuth,
-                                               dni, ghi, dhi,
-                                               dni_extra=dni_extra,
-                                               airmass=airmass,
-                                               albedo=albedo,
-                                               model=model,
-                                               **kwargs)
+        return irradiance.get_total_irradiance(
+            orientation['surface_tilt'],
+            orientation['surface_azimuth'],
+            solar_zenith, solar_azimuth,
+            dni, ghi, dhi,
+            dni_extra=dni_extra,
+            airmass=airmass,
+            albedo=albedo,
+            model=model,
+            diffuse_components=diffuse_components,
+            **kwargs
+        )
 
     def get_iam(self, aoi, iam_model='physical'):
         """
-        Determine the incidence angle modifier using the method specified by
-        ``iam_model``.
+        Determine the incidence angle modifier for direct irradiance
+        using the method specified by ``iam_model``.
 
         Parameters for the selected IAM model are expected to be in
         ``Array.module_parameters``. Default parameters are available for
@@ -1213,6 +1270,68 @@ def get_iam(self, aoi, iam_model='physical'):
         else:
             raise ValueError(model + ' is not a valid IAM model')
 
+    def get_iam_diffuse(self, surface_tilt, iam_model='marion_diffuse',
+                        marion_model=None, **kwargs):
+        """
+        Determine the incidence angle modifier for various diffuse irradiance
+        components using the method specified by ``iam_model``.
+
+        Parameters for the selected IAM model are expected to be in
+        ``Array.module_parameters``. Default parameters are available for
+        the 'marion_diffuse' and 'martin_ruiz_diffuse' models.
+
+        Parameters
+        ----------
+        surface_tilt : float or Series
+            The tilt angle of the surface in degrees.
+        iam_model : string, default 'marion_diffuse'
+            The IAM model to be used. Valid strings are 'marion_diffuse',
+            'martin_ruiz_diffuse' and 'schlick_diffuse'.
+        marion_model : string, default None
+            The IAM function to evaluate across a solid angle. Only used when
+            ``iam_model='marion_diffuse'``. Must be one of 'ashrae',
+            'physical', 'martin_ruiz' and 'sapm'.
+
+        kwargs : dict, optional
+            Additional keyword arguments passed to the IAM model function.
+
+        Returns
+        -------
+        iam_diffuse : dict or DataFrame
+            The AOI modifiers for different diffuse irradiance components.
+            Included components depend on the selected ``iam_model``.
+
+        Raises
+        ------
+        ValueError
+            if `iam_model` is not a valid model name.
+        ValueError
+            if `iam_model` is 'marion_diffuse' and `marion_model` is None.
+        """
+        model = iam_model.lower()
+        if model == 'marion_diffuse' and marion_model is None:
+            raise ValueError('marion_model must be specified when '
+                             'iam_model="marion_diffuse"')
+        if model in ['marion_diffuse', 'martin_ruiz_diffuse',
+                     'schlick_diffuse']:
+            func = getattr(iam, model)  # get function at pvlib.iam
+            # get all parameters from function signature to retrieve them from
+            # module_parameters if present
+            params = set(inspect.signature(func).parameters.keys())
+            kwargs.update(_build_kwargs(params, self.module_parameters))
+            if iam_model == 'marion_diffuse':
+                iams = func(model=marion_model, surface_tilt=surface_tilt,
+                            **kwargs)
+            else:
+                iams = func(surface_tilt=surface_tilt, **kwargs)
+        else:
+            raise ValueError(model + ' is not a valid diffuse IAM model')
+
+        if isinstance(surface_tilt, pd.Series):
+            iams = pd.DataFrame(iams, index=surface_tilt.index)
+
+        return iams
+
     def get_cell_temperature(self, poa_global, temp_air, wind_speed, model,
                              effective_irradiance=None, longwave_down=None):
         """
diff --git a/tests/test_iam.py b/tests/test_iam.py
index 123548cd6e..0287fd8fd2 100644
--- a/tests/test_iam.py
+++ b/tests/test_iam.py
@@ -134,15 +134,18 @@ def test_martin_ruiz_diffuse():
 
     surface_tilt = 30.
     a_r = 0.16
-    expected = (0.9549735, 0.7944426)
+    expected_sky = 0.9549735
+    expected_ground = 0.7944426
 
     # will fail if default values change
-    iam = _iam.martin_ruiz_diffuse(surface_tilt)
-    assert_allclose(iam, expected)
+    actual_iam = _iam.martin_ruiz_diffuse(surface_tilt)
+    assert_allclose(actual_iam['sky'], expected_sky)
+    assert_allclose(actual_iam['ground'], expected_ground)
 
     # will fail if parameter names change
     iam = _iam.martin_ruiz_diffuse(surface_tilt=surface_tilt, a_r=a_r)
-    assert_allclose(iam, expected)
+    assert_allclose(iam['sky'], expected_sky)
+    assert_allclose(iam['ground'], expected_ground)
 
     a_r = 0.18
     surface_tilt = [0, 30, 90, 120, 180, np.nan, np.inf]
@@ -153,21 +156,21 @@ def test_martin_ruiz_diffuse():
 
     # check various inputs as list
     iam = _iam.martin_ruiz_diffuse(surface_tilt, a_r)
-    assert_allclose(iam[0], expected_sky, atol=1e-7, equal_nan=True)
-    assert_allclose(iam[1], expected_gnd, atol=1e-7, equal_nan=True)
+    assert_allclose(iam['sky'], expected_sky, atol=1e-7, equal_nan=True)
+    assert_allclose(iam['ground'], expected_gnd, atol=1e-7, equal_nan=True)
 
     # check various inputs as array
     iam = _iam.martin_ruiz_diffuse(np.array(surface_tilt), a_r)
-    assert_allclose(iam[0], expected_sky, atol=1e-7, equal_nan=True)
-    assert_allclose(iam[1], expected_gnd, atol=1e-7, equal_nan=True)
+    assert_allclose(iam['sky'], expected_sky, atol=1e-7, equal_nan=True)
+    assert_allclose(iam['ground'], expected_gnd, atol=1e-7, equal_nan=True)
 
     # check various inputs as Series
     surface_tilt = pd.Series(surface_tilt)
     expected_sky = pd.Series(expected_sky, name='iam_sky')
     expected_gnd = pd.Series(expected_gnd, name='iam_ground')
     iam = _iam.martin_ruiz_diffuse(surface_tilt, a_r)
-    assert_series_equal(iam[0], expected_sky)
-    assert_series_equal(iam[1], expected_gnd)
+    assert_series_equal(iam['sky'], expected_sky)
+    assert_series_equal(iam['ground'], expected_gnd)
 
 
 def test_iam_interp():
@@ -441,22 +444,21 @@ def test_schlick_diffuse():
     expected_ground = np.array([0, 0.62693858, 0.93218737, 0.95238094])
 
     # numpy arrays
-    actual_sky, actual_ground = _iam.schlick_diffuse(surface_tilt)
-    assert_allclose(expected_sky, actual_sky)
-    assert_allclose(expected_ground, actual_ground, rtol=1e-6)
+    actual_iam = _iam.schlick_diffuse(surface_tilt)
+    assert_allclose(expected_sky, actual_iam['sky'])
+    assert_allclose(expected_ground, actual_iam['ground'], rtol=1e-6)
 
     # scalars
     for i in range(len(surface_tilt)):
-        actual_sky, actual_ground = _iam.schlick_diffuse(surface_tilt[i])
-        assert_allclose(expected_sky[i], actual_sky)
-        assert_allclose(expected_ground[i], actual_ground, rtol=1e-6)
+        actual_iam = _iam.schlick_diffuse(surface_tilt[i])
+        assert_allclose(expected_sky[i], actual_iam['sky'], rtol=1e-6)
+        assert_allclose(expected_ground[i], actual_iam['ground'], rtol=1e-6)
 
     # pandas Series
     idx = pd.date_range('2019-01-01', freq='h', periods=len(surface_tilt))
-    actual_sky, actual_ground = _iam.schlick_diffuse(pd.Series(surface_tilt,
-                                                               idx))
-    assert_series_equal(pd.Series(expected_sky, idx), actual_sky)
-    assert_series_equal(pd.Series(expected_ground, idx), actual_ground,
+    actual_iam = _iam.schlick_diffuse(pd.Series(surface_tilt, idx))
+    assert_series_equal(pd.Series(expected_sky, idx), actual_iam['sky'])
+    assert_series_equal(pd.Series(expected_ground, idx), actual_iam['ground'],
                         rtol=1e-6)
 
 
diff --git a/tests/test_pvsystem.py b/tests/test_pvsystem.py
index 78aaed70c4..90f640dfd3 100644
--- a/tests/test_pvsystem.py
+++ b/tests/test_pvsystem.py
@@ -113,6 +113,65 @@ def test_PVSystem_get_iam_invalid(sapm_module_params, mocker):
         system.get_iam(45, iam_model='not_a_model')
 
 
+def test_PVSystem_get_iam_diffuse_marion(mocker):
+    model_params = {'b': 0.05}
+    m = mocker.spy(_iam, 'marion_diffuse')
+    system = pvsystem.PVSystem(module_parameters=model_params)
+    tilt = 30
+    iam = system.get_iam_diffuse(tilt, iam_model='marion_diffuse',
+                                 marion_model='ashrae', **model_params)
+    m.assert_called_with(model='ashrae', surface_tilt=tilt, **model_params)
+    assert isinstance(iam, dict)
+    assert set(iam.keys()) == {'sky', 'ground', 'horizon'}
+
+    tilt = pd.Series([30, 60])
+    iam = system.get_iam_diffuse(tilt, iam_model='marion_diffuse',
+                                 marion_model='ashrae', **model_params)
+    assert isinstance(iam, pd.DataFrame)
+
+
+@pytest.mark.parametrize('iam_model', ['martin_ruiz_diffuse',
+                                       'schlick_diffuse'])
+def test_PVSystem_get_iam_diffuse(iam_model, mocker):
+    model_params = {'a_r': 0.16} if iam_model == 'martin_ruiz_diffuse' else {}
+    m = mocker.spy(_iam, iam_model)
+    system = pvsystem.PVSystem(module_parameters=model_params)
+    tilt = 30
+    iam = system.get_iam_diffuse(tilt, iam_model=iam_model)
+    m.assert_called_with(surface_tilt=tilt, **model_params)
+    assert isinstance(iam, dict)
+
+
+def test_PVSystem_multi_array_get_iam_diffuse():
+    model_params = {'b': 0.05}
+    system = pvsystem.PVSystem(
+        arrays=[pvsystem.Array(mount=pvsystem.FixedMount(0, 180),
+                               module_parameters=model_params),
+                pvsystem.Array(mount=pvsystem.FixedMount(0, 180),
+                               module_parameters=model_params)]
+    )
+    iam = system.get_iam_diffuse((30, 60), iam_model='marion_diffuse',
+                                 marion_model='ashrae', **model_params)
+    assert len(iam) == 2
+    assert iam[0] != iam[1]
+    with pytest.raises(ValueError,
+                       match="Length mismatch for per-array parameter"):
+        system.get_iam_diffuse((30,), iam_model='marion_diffuse',
+                               marion_model='ashrae', **model_params)
+
+
+def test_PVSystem_get_iam_diffuse_invalid(sapm_module_params):
+    system = pvsystem.PVSystem(module_parameters=sapm_module_params)
+    with pytest.raises(ValueError):
+        system.get_iam_diffuse(45, iam_model='not_a_model')
+
+
+def test_PVSystem_get_iam_diffuse_marion_missing_model(sapm_module_params):
+    system = pvsystem.PVSystem(module_parameters=sapm_module_params)
+    with pytest.raises(ValueError, match="marion_model must be specified"):
+        system.get_iam_diffuse(45, iam_model='marion_diffuse')
+
+
 def test_retrieve_sam_raises_exceptions():
     """
     Raise an exception if an invalid parameter is provided to `retrieve_sam()`.