diff --git a/docs/sphinx/source/whatsnew/v0.9.2.rst b/docs/sphinx/source/whatsnew/v0.9.2.rst
index 9a1e0b6d7b..d3e249d56b 100644
--- a/docs/sphinx/source/whatsnew/v0.9.2.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.2.rst
@@ -8,6 +8,10 @@ Deprecations
 
 Enhancements
 ~~~~~~~~~~~~
+* albedo can now be provided as a column in the `weather` DataFrame input to
+  :py:method:`pvlib.modelchain.ModelChain.run_model`. (:issue:`1387`, :pull:`1478`)
+* albedo is now available as an input to :py:meth:`pvlib.pvsystem.PVSystem.get_irradiance`
+  and :py:meth:`pvlib.pvsystem.Array.get_irradiance`. (:pull:`1478`)
 * :py:func:`pvlib.iotools.read_surfrad` now also accepts remote files
   with https links in addition to files on the SURFRAD FTP server
   (:pull:`1459`)
@@ -20,14 +24,15 @@ Enhancements
 * Add support for `PEP517 <https://peps.python.org/pep-0517/>`_ & `PEP518 <https://peps.python.org/pep-0518/>`_
   with setuptools build backend. (:pull:`1495`)
 
+
 Bug fixes
 ~~~~~~~~~
 * :py:func:`pvlib.irradiance.get_total_irradiance` and
   :py:func:`pvlib.solarposition.spa_python` now raise an error instead
-  of silently ignoring unknown parameters (:pull:`1437`)
+  of silently ignoring unknown parameters. (:pull:`1437`)
 * Fix a bug in :py:func:`pvlib.solarposition.sun_rise_set_transit_ephem`
   where passing localized timezones with large UTC offsets could return
-  rise/set/transit times for the wrong day in recent versions of ``ephem``
+  rise/set/transit times for the wrong day in recent versions of ``ephem``.
   (:issue:`1449`, :pull:`1448`)
 * :py:func:`pvlib.iotools.read_tmy3` is now able to accept midnight
   timestamps as either 24:00 (which is the standard) as well as 00:00.
@@ -68,6 +73,7 @@ Contributors
 * Naman Priyadarshi (:ghuser:`Naman-Priyadarshi`)
 * Chencheng Luo (:ghuser:`roger-lcc`)
 * Prajwal Borkar (:ghuser:`PrajwalBorkar`) 
+* Cliff Hansen (:ghuser:`cwhanse`)
 * Kevin Anderson (:ghuser:`kanderso-nrel`)
 * Cliff Hansen (:ghuser:`cwhanse`)
 * Jules Chéron (:ghuser:`jules-ch`)
diff --git a/pvlib/clearsky.py b/pvlib/clearsky.py
index 311d683674..2f5cd68c62 100644
--- a/pvlib/clearsky.py
+++ b/pvlib/clearsky.py
@@ -960,8 +960,8 @@ def bird(zenith, airmass_relative, aod380, aod500, precipitable_water,
         Extraterrestrial radiation [W/m^2], defaults to 1364[W/m^2]
     asymmetry : numeric
         Asymmetry factor, defaults to 0.85
-    albedo : numeric
-        Albedo, defaults to 0.2
+    albedo : numeric, default 0.2
+        Ground surface albedo. [unitless]
 
     Returns
     -------
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 3de4d96f65..03ddd13f5a 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -304,7 +304,7 @@ def beam_component(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
 def get_total_irradiance(surface_tilt, surface_azimuth,
                          solar_zenith, solar_azimuth,
                          dni, ghi, dhi, dni_extra=None, airmass=None,
-                         albedo=.25, surface_type=None,
+                         albedo=0.25, surface_type=None,
                          model='isotropic',
                          model_perez='allsitescomposite1990'):
     r"""
@@ -344,7 +344,7 @@ def get_total_irradiance(surface_tilt, surface_azimuth,
     airmass : None or numeric, default None
         Relative airmass (not adjusted for pressure). [unitless]
     albedo : numeric, default 0.25
-        Surface albedo. [unitless]
+        Ground surface albedo. [unitless]
     surface_type : None or str, default None
         Surface type. See :py:func:`~pvlib.irradiance.get_ground_diffuse` for
         the list of accepted values.
@@ -1872,7 +1872,7 @@ def gti_dirint(poa_global, aoi, solar_zenith, solar_azimuth, times,
         applied.
 
     albedo : numeric, default 0.25
-        Surface albedo
+        Ground surface albedo. [unitless]
 
     model : String, default 'perez'
         Irradiance model.  See :py:func:`get_sky_diffuse` for allowed values.
diff --git a/pvlib/modelchain.py b/pvlib/modelchain.py
index 2798c39a68..8211981433 100644
--- a/pvlib/modelchain.py
+++ b/pvlib/modelchain.py
@@ -268,7 +268,7 @@ class ModelChainResult:
     _per_array_fields = {'total_irrad', 'aoi', 'aoi_modifier',
                          'spectral_modifier', 'cell_temperature',
                          'effective_irradiance', 'dc', 'diode_params',
-                         'dc_ohmic_losses', 'weather'}
+                         'dc_ohmic_losses', 'weather', 'albedo'}
 
     # system-level information
     solar_position: Optional[pd.DataFrame] = field(default=None)
@@ -366,6 +366,10 @@ class ModelChainResult:
     """DatetimeIndex containing a copy of the index of the input weather data.
     """
 
+    albedo: Optional[PerArray[pd.Series]] = None
+    """Series (or tuple of Series, one for each array) containing albedo.
+    """
+
     def _result_type(self, value):
         """Coerce `value` to the correct type according to
         ``self._singleton_tuples``."""
@@ -1339,6 +1343,17 @@ def _prep_inputs_solar_pos(self, weather):
             **kwargs)
         return self
 
+    def _prep_inputs_albedo(self, weather):
+        """
+        Get albedo from weather
+        """
+        try:
+            self.results.albedo = _tuple_from_dfs(weather, 'albedo')
+        except KeyError:
+            self.results.albedo = tuple([
+                a.albedo for a in self.system.arrays])
+        return self
+
     def _prep_inputs_airmass(self):
         """
         Assign airmass
@@ -1471,11 +1486,17 @@ def prepare_inputs(self, weather):
 
         Parameters
         ----------
-        weather : DataFrame, or tuple or list of DataFrame
+        weather : DataFrame, or tuple or list of DataFrames
             Required column names include ``'dni'``, ``'ghi'``, ``'dhi'``.
-            Optional column names are ``'wind_speed'``, ``'temp_air'``; if not
+            Optional column names are ``'wind_speed'``, ``'temp_air'``,
+            ``'albedo'``.
+
+            If optional columns ``'wind_speed'``, ``'temp_air'`` are not
             provided, air temperature of 20 C and wind speed
-            of 0 m/s will be added to the DataFrame.
+            of 0 m/s will be added to the ``weather`` DataFrame.
+
+            If optional column ``'albedo'`` is provided, albedo values in the
+            ModelChain's PVSystem.arrays are ignored.
 
             If `weather` is a tuple or list, it must be of the same length and
             order as the Arrays of the ModelChain's PVSystem.
@@ -1494,7 +1515,7 @@ def prepare_inputs(self, weather):
         Notes
         -----
         Assigns attributes to ``results``: ``times``, ``weather``,
-        ``solar_position``, ``airmass``, ``total_irrad``, ``aoi``
+        ``solar_position``, ``airmass``, ``total_irrad``, ``aoi``, ``albedo``.
 
         See also
         --------
@@ -1507,6 +1528,7 @@ def prepare_inputs(self, weather):
 
         self._prep_inputs_solar_pos(weather)
         self._prep_inputs_airmass()
+        self._prep_inputs_albedo(weather)
 
         # PVSystem.get_irradiance and SingleAxisTracker.get_irradiance
         # and PVSystem.get_aoi and SingleAxisTracker.get_aoi
@@ -1531,6 +1553,7 @@ def prepare_inputs(self, weather):
             _tuple_from_dfs(self.results.weather, 'dni'),
             _tuple_from_dfs(self.results.weather, 'ghi'),
             _tuple_from_dfs(self.results.weather, 'dhi'),
+            albedo=self.results.albedo,
             airmass=self.results.airmass['airmass_relative'],
             model=self.transposition_model
         )
@@ -1724,16 +1747,32 @@ def run_model(self, weather):
         Parameters
         ----------
         weather : DataFrame, or tuple or list of DataFrame
-            Irradiance column names must include ``'dni'``, ``'ghi'``, and
-            ``'dhi'``. If optional columns ``'temp_air'`` and ``'wind_speed'``
+            Column names must include:
+
+            - ``'dni'``
+            - ``'ghi'``
+            - ``'dhi'``
+
+            Optional columns are:
+
+            - ``'temp_air'``
+            - ``'cell_temperature'``
+            - ``'module_temperature'``
+            - ``'wind_speed'``
+            - ``'albedo'``
+
+            If optional columns ``'temp_air'`` and ``'wind_speed'``
             are not provided, air temperature of 20 C and wind speed of 0 m/s
             are added to the DataFrame. If optional column
             ``'cell_temperature'`` is provided, these values are used instead
-            of `temperature_model`. If optional column `module_temperature`
-            is provided, `temperature_model` must be ``'sapm'``.
+            of `temperature_model`. If optional column ``'module_temperature'``
+            is provided, ``temperature_model`` must be ``'sapm'``.
 
-            If list or tuple, must be of the same length and order as the
-            Arrays of the ModelChain's PVSystem.
+            If optional column ``'albedo'`` is provided, ``'albedo'`` may not
+            be present on the ModelChain's PVSystem.Arrays.
+
+            If weather is a list or tuple, it must be of the same length and
+            order as the Arrays of the ModelChain's PVSystem.
 
         Returns
         -------
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 6bb89f34a3..77560e04c0 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -134,7 +134,7 @@ class PVSystem:
         a single array is created from the other parameters (e.g.
         `surface_tilt`, `surface_azimuth`). Must contain at least one Array,
         if length of arrays is 0 a ValueError is raised. If `arrays` is
-        specified the following parameters are ignored:
+        specified the following PVSystem parameters are ignored:
 
         - `surface_tilt`
         - `surface_azimuth`
@@ -157,13 +157,14 @@ class PVSystem:
         North=0, East=90, South=180, West=270.
 
     albedo : None or float, default None
-        The ground albedo. If ``None``, will attempt to use
-        ``surface_type`` and ``irradiance.SURFACE_ALBEDOS``
-        to lookup albedo.
+        Ground surface albedo. If ``None``, then ``surface_type`` is used
+        to look up a value in ``irradiance.SURFACE_ALBEDOS``.
+        If ``surface_type`` is also None then a ground surface albedo
+        of 0.25 is used.
 
     surface_type : None or string, default None
-        The ground surface type. See ``irradiance.SURFACE_ALBEDOS``
-        for valid values.
+        The ground surface type. See ``irradiance.SURFACE_ALBEDOS`` for
+        valid values.
 
     module : None or string, default None
         The model name of the modules.
@@ -333,30 +334,33 @@ 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, model='haydavies',
-                       **kwargs):
+                       dni_extra=None, airmass=None, albedo=None,
+                       model='haydavies', **kwargs):
         """
         Uses the :py:func:`irradiance.get_total_irradiance` function to
-        calculate the plane of array irradiance components on a tilted
-        surface defined by ``self.surface_tilt``,
-        ``self.surface_azimuth``, and ``self.albedo``.
+        calculate the plane of array irradiance components on the tilted
+        surfaces defined by each array's ``surface_tilt`` and
+        ``surface_azimuth``.
 
         Parameters
         ----------
-        solar_zenith : float or Series.
+        solar_zenith : float or Series
             Solar zenith angle.
-        solar_azimuth : float or Series.
+        solar_azimuth : float or Series
             Solar azimuth angle.
         dni : float or Series or tuple of float or Series
-            Direct Normal Irradiance
+            Direct Normal Irradiance. [W/m2]
         ghi : float or Series or tuple of float or Series
-            Global horizontal irradiance
+            Global horizontal irradiance. [W/m2]
         dhi : float or Series or tuple of float or Series
-            Diffuse horizontal irradiance
-        dni_extra : None, float or Series, default None
-            Extraterrestrial direct normal irradiance
+            Diffuse horizontal irradiance. [W/m2]
+        dni_extra : None, float, Series or tuple of float or Series,\
+            default None
+            Extraterrestrial direct normal irradiance. [W/m2]
         airmass : None, float or Series, default None
-            Airmass
+            Airmass. [unitless]
+        albedo : None, float or Series, default None
+            Ground surface albedo. [unitless]
         model : String, default 'haydavies'
             Irradiance model.
 
@@ -376,17 +380,24 @@ def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
         poa_irradiance : DataFrame or tuple of DataFrame
             Column names are: ``'poa_global', 'poa_direct', 'poa_diffuse',
             'poa_sky_diffuse', 'poa_ground_diffuse'``.
+
+        See also
+        --------
+        pvlib.irradiance.get_total_irradiance
         """
         dni = self._validate_per_array(dni, system_wide=True)
         ghi = self._validate_per_array(ghi, system_wide=True)
         dhi = self._validate_per_array(dhi, system_wide=True)
+
+        albedo = self._validate_per_array(albedo, system_wide=True)
+
         return tuple(
             array.get_irradiance(solar_zenith, solar_azimuth,
                                  dni, ghi, dhi,
-                                 dni_extra, airmass, model,
-                                 **kwargs)
-            for array, dni, ghi, dhi in zip(
-                self.arrays, dni, ghi, dhi
+                                 dni_extra=dni_extra, airmass=airmass,
+                                 albedo=albedo, model=model, **kwargs)
+            for array, dni, ghi, dhi, albedo in zip(
+                self.arrays, dni, ghi, dhi, albedo
             )
         )
 
@@ -1258,14 +1269,14 @@ class Array:
         If not provided, a FixedMount with zero tilt is used.
 
     albedo : None or float, default None
-        The ground albedo. If ``None``, will attempt to use
-        ``surface_type`` to look up an albedo value in
-        ``irradiance.SURFACE_ALBEDOS``. If a surface albedo
-        cannot be found then 0.25 is used.
+        Ground surface albedo. If ``None``, then ``surface_type`` is used
+        to look up a value in ``irradiance.SURFACE_ALBEDOS``.
+        If ``surface_type`` is also None then a ground surface albedo
+        of 0.25 is used.
 
     surface_type : None or string, default None
-        The ground surface type. See ``irradiance.SURFACE_ALBEDOS``
-        for valid values.
+        The ground surface type. See ``irradiance.SURFACE_ALBEDOS`` for valid
+        values.
 
     module : None or string, default None
         The model name of the modules.
@@ -1425,15 +1436,14 @@ def get_aoi(self, solar_zenith, solar_azimuth):
                               solar_zenith, solar_azimuth)
 
     def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
-                       dni_extra=None, airmass=None, model='haydavies',
-                       **kwargs):
+                       dni_extra=None, airmass=None, albedo=None,
+                       model='haydavies', **kwargs):
         """
         Get plane of array irradiance components.
 
         Uses the :py:func:`pvlib.irradiance.get_total_irradiance` function to
         calculate the plane of array irradiance components for a surface
-        defined by ``self.surface_tilt`` and ``self.surface_azimuth`` with
-        albedo ``self.albedo``.
+        defined by ``self.surface_tilt`` and ``self.surface_azimuth``.
 
         Parameters
         ----------
@@ -1442,15 +1452,17 @@ def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
         solar_azimuth : float or Series.
             Solar azimuth angle.
         dni : float or Series
-            Direct Normal Irradiance
-        ghi : float or Series
+            Direct normal irradiance. [W/m2]
+        ghi : float or Series. [W/m2]
             Global horizontal irradiance
         dhi : float or Series
-            Diffuse horizontal irradiance
+            Diffuse horizontal irradiance. [W/m2]
         dni_extra : None, float or Series, default None
-            Extraterrestrial direct normal irradiance
+            Extraterrestrial direct normal irradiance. [W/m2]
         airmass : None, float or Series, default None
-            Airmass
+            Airmass. [unitless]
+        albedo : None, float or Series, default None
+            Ground surface albedo. [unitless]
         model : String, default 'haydavies'
             Irradiance model.
 
@@ -1463,7 +1475,14 @@ def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
         poa_irradiance : DataFrame
             Column names are: ``'poa_global', 'poa_direct', 'poa_diffuse',
             'poa_sky_diffuse', 'poa_ground_diffuse'``.
+
+        See also
+        --------
+        :py:func:`pvlib.irradiance.get_total_irradiance`
         """
+        if albedo is None:
+            albedo = self.albedo
+
         # not needed for all models, but this is easier
         if dni_extra is None:
             dni_extra = irradiance.get_extra_radiation(solar_zenith.index)
@@ -1478,8 +1497,8 @@ def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
                                                dni, ghi, dhi,
                                                dni_extra=dni_extra,
                                                airmass=airmass,
+                                               albedo=albedo,
                                                model=model,
-                                               albedo=self.albedo,
                                                **kwargs)
 
     def get_iam(self, aoi, iam_model='physical'):
@@ -3293,7 +3312,7 @@ def dc_ohms_from_percent(vmp_ref, imp_ref, dc_ohmic_percent,
 
     See Also
     --------
-    :py:func:`~pvlib.pvsystem.dc_ohmic_losses`
+    pvlib.pvsystem.dc_ohmic_losses
 
     References
     ----------
@@ -3328,7 +3347,7 @@ def dc_ohmic_losses(resistance, current):
 
     See Also
     --------
-    :py:func:`~pvlib.pvsystem.dc_ohms_from_percent`
+    pvlib.pvsystem.dc_ohms_from_percent
 
     References
     ----------
diff --git a/pvlib/tests/test_clearsky.py b/pvlib/tests/test_clearsky.py
index 15fc74e383..d603cbcdfe 100644
--- a/pvlib/tests/test_clearsky.py
+++ b/pvlib/tests/test_clearsky.py
@@ -756,6 +756,30 @@ def test_bird():
     assert np.allclose(
         testdata['Dif Hz'].where(dusk, 0.), diffuse_horz[1:48], rtol=1e-3
     )
+    # repeat test with albedo as a Series
+    alb_series = pd.Series(0.2, index=times)
+    irrads = clearsky.bird(
+        zenith, airmass, aod_380nm, aod_500nm, h2o_cm, o3_cm, press_mB * 100.,
+        etr, b_a, alb_series
+    )
+    Eb, Ebh, Gh, Dh = (irrads[_] for _ in field_names)
+    direct_beam = pd.Series(np.where(dawn, Eb, 0.), index=times).fillna(0.)
+    assert np.allclose(
+        testdata['Direct Beam'].where(dusk, 0.), direct_beam[1:48], rtol=1e-3
+    )
+    direct_horz = pd.Series(np.where(dawn, Ebh, 0.), index=times).fillna(0.)
+    assert np.allclose(
+        testdata['Direct Hz'].where(dusk, 0.), direct_horz[1:48], rtol=1e-3
+    )
+    global_horz = pd.Series(np.where(dawn, Gh, 0.), index=times).fillna(0.)
+    assert np.allclose(
+        testdata['Global Hz'].where(dusk, 0.), global_horz[1:48], rtol=1e-3
+    )
+    diffuse_horz = pd.Series(np.where(dawn, Dh, 0.), index=times).fillna(0.)
+    assert np.allclose(
+        testdata['Dif Hz'].where(dusk, 0.), diffuse_horz[1:48], rtol=1e-3
+    )
+
     # test keyword parameters
     irrads2 = clearsky.bird(
         zenith, airmass, aod_380nm, aod_500nm, h2o_cm, dni_extra=etr
diff --git a/pvlib/tests/test_irradiance.py b/pvlib/tests/test_irradiance.py
index 80986f26c3..8dc4877d0d 100644
--- a/pvlib/tests/test_irradiance.py
+++ b/pvlib/tests/test_irradiance.py
@@ -120,29 +120,38 @@ def test_get_extra_radiation_invalid():
         irradiance.get_extra_radiation(300, method='invalid')
 
 
-def test_grounddiffuse_simple_float():
+def test_get_ground_diffuse_simple_float():
     result = irradiance.get_ground_diffuse(40, 900)
     assert_allclose(result, 26.32000014911496)
 
 
-def test_grounddiffuse_simple_series(irrad_data):
+def test_get_ground_diffuse_simple_series(irrad_data):
     ground_irrad = irradiance.get_ground_diffuse(40, irrad_data['ghi'])
     assert ground_irrad.name == 'diffuse_ground'
 
 
-def test_grounddiffuse_albedo_0(irrad_data):
+def test_get_ground_diffuse_albedo_0(irrad_data):
     ground_irrad = irradiance.get_ground_diffuse(
         40, irrad_data['ghi'], albedo=0)
     assert 0 == ground_irrad.all()
 
 
+def test_get_ground_diffuse_albedo_series(times):
+    albedo = pd.Series(0.2, index=times)
+    ground_irrad = irradiance.get_ground_diffuse(
+        45, pd.Series(1000, index=times), albedo)
+    expected = albedo * 0.5 * (1 - np.sqrt(2) / 2.) * 1000
+    expected.name = 'diffuse_ground'
+    assert_series_equal(ground_irrad, expected)
+
+
 def test_grounddiffuse_albedo_invalid_surface(irrad_data):
     with pytest.raises(KeyError):
         irradiance.get_ground_diffuse(
             40, irrad_data['ghi'], surface_type='invalid')
 
 
-def test_grounddiffuse_albedo_surface(irrad_data):
+def test_get_ground_diffuse_albedo_surface(irrad_data):
     result = irradiance.get_ground_diffuse(40, irrad_data['ghi'],
                                            surface_type='sand')
     assert_allclose(result, [0, 3.731058, 48.778813, 12.035025], atol=1e-4)
@@ -387,6 +396,25 @@ def test_get_total_irradiance(irrad_data, ephem_data, dni_et,
                                           'poa_ground_diffuse']
 
 
+@pytest.mark.parametrize('model', ['isotropic', 'klucher',
+                                   'haydavies', 'reindl', 'king', 'perez'])
+def test_get_total_irradiance_albedo(
+        irrad_data, ephem_data, dni_et, relative_airmass, model):
+    albedo = pd.Series(0.2, index=ephem_data.index)
+    total = irradiance.get_total_irradiance(
+        32, 180,
+        ephem_data['apparent_zenith'], ephem_data['azimuth'],
+        dni=irrad_data['dni'], ghi=irrad_data['ghi'],
+        dhi=irrad_data['dhi'],
+        dni_extra=dni_et, airmass=relative_airmass,
+        model=model,
+        albedo=albedo)
+
+    assert total.columns.tolist() == ['poa_global', 'poa_direct',
+                                      'poa_diffuse', 'poa_sky_diffuse',
+                                      'poa_ground_diffuse']
+
+
 @pytest.mark.parametrize('model', ['isotropic', 'klucher',
                                    'haydavies', 'reindl', 'king', 'perez'])
 def test_get_total_irradiance_scalars(model):
@@ -698,6 +726,14 @@ def test_gti_dirint():
 
     assert_frame_equal(output, expected)
 
+    # test with albedo as a Series
+    albedo = pd.Series(0.05, index=times)
+    output = irradiance.gti_dirint(
+        poa_global, aoi, zenith, azimuth, times, surface_tilt, surface_azimuth,
+        albedo=albedo)
+
+    assert_frame_equal(output, expected)
+
     # test temp_dew input
     temp_dew = np.array([70, 80, 20])
     output = irradiance.gti_dirint(
diff --git a/pvlib/tests/test_modelchain.py b/pvlib/tests/test_modelchain.py
index f4a92eadad..62b71f2042 100644
--- a/pvlib/tests/test_modelchain.py
+++ b/pvlib/tests/test_modelchain.py
@@ -495,6 +495,26 @@ def test_prepare_inputs_multi_weather(
     mc.prepare_inputs(input_type((weather, weather)))
     num_arrays = sapm_dc_snl_ac_system_Array.num_arrays
     assert len(mc.results.total_irrad) == num_arrays
+    # check that albedo is transfered to mc.results from mc.system.arrays
+    assert mc.results.albedo == (0.2, 0.2)
+
+
+@pytest.mark.parametrize("input_type", [tuple, list])
+def test_prepare_inputs_albedo_in_weather(
+        sapm_dc_snl_ac_system_Array, location, input_type):
+    times = pd.date_range(start='20160101 1200-0700',
+                          end='20160101 1800-0700', freq='6H')
+    mc = ModelChain(sapm_dc_snl_ac_system_Array, location)
+    weather = pd.DataFrame({'ghi': 1, 'dhi': 1, 'dni': 1, 'albedo': 0.5},
+                           index=times)
+    # weather as a single DataFrame
+    mc.prepare_inputs(weather)
+    num_arrays = sapm_dc_snl_ac_system_Array.num_arrays
+    assert len(mc.results.albedo) == num_arrays
+    # repeat with tuple of weather
+    mc.prepare_inputs(input_type((weather, weather)))
+    num_arrays = sapm_dc_snl_ac_system_Array.num_arrays
+    assert len(mc.results.albedo) == num_arrays
 
 
 def test_prepare_inputs_no_irradiance(sapm_dc_snl_ac_system, location):
diff --git a/pvlib/tests/test_pvsystem.py b/pvlib/tests/test_pvsystem.py
index 1141e490e9..7fa013d0dc 100644
--- a/pvlib/tests/test_pvsystem.py
+++ b/pvlib/tests/test_pvsystem.py
@@ -18,6 +18,7 @@
 from pvlib.pvsystem import FixedMount
 from pvlib import temperature
 from pvlib._deprecation import pvlibDeprecationWarning
+from pvlib.tools import cosd
 
 
 @pytest.mark.parametrize('iam_model,model_params', [
@@ -1673,51 +1674,70 @@ def test_PVSystem_multiple_array_get_aoi():
     assert aoi_one > 0
 
 
-def test_PVSystem_get_irradiance():
-    system = pvsystem.PVSystem(surface_tilt=32, surface_azimuth=135)
+@pytest.fixture
+def solar_pos():
     times = pd.date_range(start='20160101 1200-0700',
                           end='20160101 1800-0700', freq='6H')
     location = Location(latitude=32, longitude=-111)
-    solar_position = location.get_solarposition(times)
+    return location.get_solarposition(times)
+
+
+def test_PVSystem_get_irradiance(solar_pos):
+    system = pvsystem.PVSystem(surface_tilt=32, surface_azimuth=135)
     irrads = pd.DataFrame({'dni':[900,0], 'ghi':[600,0], 'dhi':[100,0]},
-                          index=times)
+                          index=solar_pos.index)
 
-    irradiance = system.get_irradiance(solar_position['apparent_zenith'],
-                                       solar_position['azimuth'],
+    irradiance = system.get_irradiance(solar_pos['apparent_zenith'],
+                                       solar_pos['azimuth'],
                                        irrads['dni'],
                                        irrads['ghi'],
                                        irrads['dhi'])
 
     expected = pd.DataFrame(data=np.array(
-        [[ 883.65494055,  745.86141676,  137.79352379,  126.397131  ,
-              11.39639279],
-           [   0.        ,   -0.        ,    0.        ,    0.        ,    0.        ]]),
+        [[883.65494055, 745.86141676, 137.79352379, 126.397131, 11.39639279],
+         [0., -0., 0., 0., 0.]]),
                             columns=['poa_global', 'poa_direct',
                                      'poa_diffuse', 'poa_sky_diffuse',
                                      'poa_ground_diffuse'],
-                            index=times)
+                            index=solar_pos.index)
+    assert_frame_equal(irradiance, expected, check_less_precise=2)
+
 
+def test_PVSystem_get_irradiance_albedo(solar_pos):
+    system = pvsystem.PVSystem(surface_tilt=32, surface_azimuth=135)
+    irrads = pd.DataFrame({'dni': [900, 0], 'ghi': [600, 0], 'dhi': [100, 0],
+                           'albedo': [0.5, 0.5]},
+                          index=solar_pos.index)
+    # albedo as a Series
+    irradiance = system.get_irradiance(solar_pos['apparent_zenith'],
+                                       solar_pos['azimuth'],
+                                       irrads['dni'],
+                                       irrads['ghi'],
+                                       irrads['dhi'],
+                                       albedo=irrads['albedo'])
+    expected = pd.DataFrame(data=np.array(
+        [[895.05134334, 745.86141676, 149.18992658, 126.397131, 22.79279558],
+         [0., -0., 0., 0., 0.]]),
+        columns=['poa_global', 'poa_direct', 'poa_diffuse', 'poa_sky_diffuse',
+                 'poa_ground_diffuse'],
+        index=solar_pos.index)
     assert_frame_equal(irradiance, expected, check_less_precise=2)
 
 
-def test_PVSystem_get_irradiance_model(mocker):
+def test_PVSystem_get_irradiance_model(mocker, solar_pos):
     spy_perez = mocker.spy(irradiance, 'perez')
     spy_haydavies = mocker.spy(irradiance, 'haydavies')
     system = pvsystem.PVSystem(surface_tilt=32, surface_azimuth=135)
-    times = pd.date_range(start='20160101 1200-0700',
-                          end='20160101 1800-0700', freq='6H')
-    location = Location(latitude=32, longitude=-111)
-    solar_position = location.get_solarposition(times)
     irrads = pd.DataFrame({'dni': [900, 0], 'ghi': [600, 0], 'dhi': [100, 0]},
-                          index=times)
-    system.get_irradiance(solar_position['apparent_zenith'],
-                          solar_position['azimuth'],
+                          index=solar_pos.index)
+    system.get_irradiance(solar_pos['apparent_zenith'],
+                          solar_pos['azimuth'],
                           irrads['dni'],
                           irrads['ghi'],
                           irrads['dhi'])
     spy_haydavies.assert_called_once()
-    system.get_irradiance(solar_position['apparent_zenith'],
-                          solar_position['azimuth'],
+    system.get_irradiance(solar_pos['apparent_zenith'],
+                          solar_pos['azimuth'],
                           irrads['dni'],
                           irrads['ghi'],
                           irrads['dhi'],
@@ -1725,31 +1745,28 @@ def test_PVSystem_get_irradiance_model(mocker):
     spy_perez.assert_called_once()
 
 
-def test_PVSystem_multi_array_get_irradiance():
+def test_PVSystem_multi_array_get_irradiance(solar_pos):
     array_one = pvsystem.Array(pvsystem.FixedMount(surface_tilt=32,
                                                    surface_azimuth=135))
     array_two = pvsystem.Array(pvsystem.FixedMount(surface_tilt=5,
                                                    surface_azimuth=150))
     system = pvsystem.PVSystem(arrays=[array_one, array_two])
-    location = Location(latitude=32, longitude=-111)
-    times = pd.date_range(start='20160101 1200-0700',
-                          end='20160101 1800-0700', freq='6H')
-    solar_position = location.get_solarposition(times)
+
     irrads = pd.DataFrame({'dni': [900, 0], 'ghi': [600, 0], 'dhi': [100, 0]},
-                          index=times)
+                          index=solar_pos.index)
     array_one_expected = array_one.get_irradiance(
-        solar_position['apparent_zenith'],
-        solar_position['azimuth'],
+        solar_pos['apparent_zenith'],
+        solar_pos['azimuth'],
         irrads['dni'], irrads['ghi'], irrads['dhi']
     )
     array_two_expected = array_two.get_irradiance(
-        solar_position['apparent_zenith'],
-        solar_position['azimuth'],
+        solar_pos['apparent_zenith'],
+        solar_pos['azimuth'],
         irrads['dni'], irrads['ghi'], irrads['dhi']
     )
     array_one_irrad, array_two_irrad = system.get_irradiance(
-        solar_position['apparent_zenith'],
-        solar_position['azimuth'],
+        solar_pos['apparent_zenith'],
+        solar_pos['azimuth'],
         irrads['dni'], irrads['ghi'], irrads['dhi']
     )
     assert_frame_equal(
@@ -1760,7 +1777,7 @@ def test_PVSystem_multi_array_get_irradiance():
     )
 
 
-def test_PVSystem_multi_array_get_irradiance_multi_irrad():
+def test_PVSystem_multi_array_get_irradiance_multi_irrad(solar_pos):
     """Test a system with two identical arrays but different irradiance.
 
     Because only the irradiance is different we expect the same output
@@ -1771,39 +1788,36 @@ def test_PVSystem_multi_array_get_irradiance_multi_irrad():
     array_one = pvsystem.Array(pvsystem.FixedMount(0, 180))
     array_two = pvsystem.Array(pvsystem.FixedMount(0, 180))
     system = pvsystem.PVSystem(arrays=[array_one, array_two])
-    location = Location(latitude=32, longitude=-111)
-    times = pd.date_range(start='20160101 1200-0700',
-                          end='20160101 1800-0700', freq='6H')
-    solar_position = location.get_solarposition(times)
+
     irrads = pd.DataFrame({'dni': [900, 0], 'ghi': [600, 0], 'dhi': [100, 0]},
-                          index=times)
+                          index=solar_pos.index)
     irrads_two = pd.DataFrame(
         {'dni': [0, 900], 'ghi': [0, 600], 'dhi': [0, 100]},
-        index=times
+        index=solar_pos.index
     )
     array_irrad = system.get_irradiance(
-        solar_position['apparent_zenith'],
-        solar_position['azimuth'],
+        solar_pos['apparent_zenith'],
+        solar_pos['azimuth'],
         (irrads['dhi'], irrads['dhi']),
         (irrads['ghi'], irrads['ghi']),
         (irrads['dni'], irrads['dni'])
     )
     assert_frame_equal(array_irrad[0], array_irrad[1])
     array_irrad = system.get_irradiance(
-        solar_position['apparent_zenith'],
-        solar_position['azimuth'],
+        solar_pos['apparent_zenith'],
+        solar_pos['azimuth'],
         (irrads['dhi'], irrads_two['dhi']),
         (irrads['ghi'], irrads_two['ghi']),
         (irrads['dni'], irrads_two['dni'])
     )
     array_one_expected = array_one.get_irradiance(
-        solar_position['apparent_zenith'],
-        solar_position['azimuth'],
+        solar_pos['apparent_zenith'],
+        solar_pos['azimuth'],
         irrads['dhi'], irrads['ghi'], irrads['dni']
     )
     array_two_expected = array_two.get_irradiance(
-        solar_position['apparent_zenith'],
-        solar_position['azimuth'],
+        solar_pos['apparent_zenith'],
+        solar_pos['azimuth'],
         irrads_two['dhi'], irrads_two['ghi'], irrads_two['dni']
     )
     assert not array_irrad[0].equals(array_irrad[1])
@@ -1812,15 +1826,15 @@ def test_PVSystem_multi_array_get_irradiance_multi_irrad():
     with pytest.raises(ValueError,
                        match="Length mismatch for per-array parameter"):
         system.get_irradiance(
-            solar_position['apparent_zenith'],
-            solar_position['azimuth'],
+            solar_pos['apparent_zenith'],
+            solar_pos['azimuth'],
             (irrads['dhi'], irrads_two['dhi'], irrads['dhi']),
             (irrads['ghi'], irrads_two['ghi']),
             irrads['dni']
         )
     array_irrad = system.get_irradiance(
-        solar_position['apparent_zenith'],
-        solar_position['azimuth'],
+        solar_pos['apparent_zenith'],
+        solar_pos['azimuth'],
         (irrads['dhi'], irrads_two['dhi']),
         irrads['ghi'],
         irrads['dni']
@@ -1829,6 +1843,44 @@ def test_PVSystem_multi_array_get_irradiance_multi_irrad():
     assert not array_irrad[0].equals(array_irrad[1])
 
 
+def test_Array_get_irradiance(solar_pos):
+    array = pvsystem.Array(pvsystem.FixedMount(surface_tilt=32,
+                                               surface_azimuth=135))
+    irrads = pd.DataFrame({'dni': [900, 0], 'ghi': [600, 0], 'dhi': [100, 0]},
+                          index=solar_pos.index)
+    # defaults for kwargs
+    modeled = array.get_irradiance(
+        solar_pos['apparent_zenith'],
+        solar_pos['azimuth'],
+        irrads['dni'], irrads['ghi'], irrads['dhi']
+    )
+    expected = pd.DataFrame(
+        data=np.array(
+            [[883.65494055, 745.86141676, 137.79352379, 126.397131,
+              11.39639279],
+             [0., -0., 0., 0., 0.]]),
+        columns=['poa_global', 'poa_direct', 'poa_diffuse', 'poa_sky_diffuse',
+                 'poa_ground_diffuse'],
+        index=solar_pos.index
+    )
+    assert_frame_equal(modeled, expected, check_less_precise=5)
+    # with specified kwargs, use isotropic sky diffuse because it's easier
+    modeled = array.get_irradiance(
+        solar_pos['apparent_zenith'],
+        solar_pos['azimuth'],
+        irrads['dni'], irrads['ghi'], irrads['dhi'],
+        albedo=0.5, model='isotropic'
+    )
+    sky_diffuse = irradiance.isotropic(array.mount.surface_tilt, irrads['dhi'])
+    ground_diff = irradiance.get_ground_diffuse(
+        array.mount.surface_tilt, irrads['ghi'], 0.5, surface_type=None)
+    aoi = irradiance.aoi(array.mount.surface_tilt, array.mount.surface_azimuth,
+                         solar_pos['apparent_zenith'], solar_pos['azimuth'])
+    direct = irrads['dni'] * cosd(aoi)
+    expected = sky_diffuse + ground_diff + direct
+    assert_series_equal(expected, expected, check_less_precise=5)
+
+
 @fail_on_pvlib_version('0.10')
 @pytest.mark.parametrize('attr', ['module_parameters', 'module', 'module_type',
                                   'temperature_model_parameters', 'albedo',
diff --git a/pvlib/tests/test_tracking.py b/pvlib/tests/test_tracking.py
index 0fbace0f17..87452939f5 100644
--- a/pvlib/tests/test_tracking.py
+++ b/pvlib/tests/test_tracking.py
@@ -393,6 +393,25 @@ def test_get_irradiance():
 
     assert_frame_equal(irradiance, expected, check_less_precise=2)
 
+    # test with albedo as a Series
+    irrads['albedo'] = [0.5, 0.5]
+    with np.errstate(invalid='ignore'):
+        irradiance = system.get_irradiance(tracker_data['surface_tilt'],
+                                           tracker_data['surface_azimuth'],
+                                           solar_zenith,
+                                           solar_azimuth,
+                                           irrads['dni'],
+                                           irrads['ghi'],
+                                           irrads['dhi'],
+                                           albedo=irrads['albedo'])
+
+    expected = pd.Series(data=[21.05514984, nan], index=times,
+                         name='poa_ground_diffuse')
+
+    assert_series_equal(irradiance['poa_ground_diffuse'], expected,
+                        check_less_precise=2)
+
+
 
 def test_SingleAxisTracker___repr__():
     with pytest.warns(pvlibDeprecationWarning):
diff --git a/pvlib/tracking.py b/pvlib/tracking.py
index 9bd216207d..7afb46f253 100644
--- a/pvlib/tracking.py
+++ b/pvlib/tracking.py
@@ -188,7 +188,8 @@ def get_aoi(self, surface_tilt, surface_azimuth, solar_zenith,
     @_unwrap_single_value
     def get_irradiance(self, surface_tilt, surface_azimuth,
                        solar_zenith, solar_azimuth, dni, ghi, dhi,
-                       dni_extra=None, airmass=None, model='haydavies',
+                       albedo=None, dni_extra=None, airmass=None,
+                       model='haydavies',
                        **kwargs):
         """
         Uses the :func:`irradiance.get_total_irradiance` function to
@@ -215,6 +216,8 @@ def get_irradiance(self, surface_tilt, surface_azimuth,
             Global horizontal irradiance
         dhi : float or Series
             Diffuse horizontal irradiance
+        albedo : None, float or Series, default None
+            Ground surface albedo. [unitless]
         dni_extra : float or Series, default None
             Extraterrestrial direct normal irradiance
         airmass : float or Series, default None
@@ -245,6 +248,13 @@ def get_irradiance(self, surface_tilt, surface_azimuth,
         ghi = self._validate_per_array(ghi, system_wide=True)
         dhi = self._validate_per_array(dhi, system_wide=True)
 
+        if albedo is None:
+            # assign default albedo here because SingleAxisTracker
+            # initializes albedo to None
+            albedo = 0.25
+
+        albedo = self._validate_per_array(albedo, system_wide=True)
+
         return tuple(
             irradiance.get_total_irradiance(
                 surface_tilt,
@@ -255,10 +265,10 @@ def get_irradiance(self, surface_tilt, surface_azimuth,
                 dni_extra=dni_extra,
                 airmass=airmass,
                 model=model,
-                albedo=self.arrays[0].albedo,
+                albedo=albedo,
                 **kwargs)
-            for array, dni, ghi, dhi in zip(
-                self.arrays, dni, ghi, dhi
+            for array, dni, ghi, dhi, albedo in zip(
+                self.arrays, dni, ghi, dhi, albedo
             )
         )