diff --git a/docs/sphinx/source/reference/pv_modeling.rst b/docs/sphinx/source/reference/pv_modeling.rst
index 31c380c1bb..0f33cf8c70 100644
--- a/docs/sphinx/source/reference/pv_modeling.rst
+++ b/docs/sphinx/source/reference/pv_modeling.rst
@@ -28,6 +28,8 @@ Incident angle modifiers
    iam.interp
    iam.marion_diffuse
    iam.marion_integrate
+   iam.schlick
+   iam.schlick_diffuse
 
 PV temperature models
 ---------------------
diff --git a/docs/sphinx/source/whatsnew.rst b/docs/sphinx/source/whatsnew.rst
index 4830371985..464e59f121 100644
--- a/docs/sphinx/source/whatsnew.rst
+++ b/docs/sphinx/source/whatsnew.rst
@@ -6,6 +6,7 @@ What's New
 
 These are new features and improvements of note in each release.
 
+.. include:: whatsnew/v0.9.4.rst
 .. include:: whatsnew/v0.9.3.rst
 .. include:: whatsnew/v0.9.2.rst
 .. include:: whatsnew/v0.9.1.rst
diff --git a/docs/sphinx/source/whatsnew/v0.9.4.rst b/docs/sphinx/source/whatsnew/v0.9.4.rst
index 8a67c201f0..93d056e5a7 100644
--- a/docs/sphinx/source/whatsnew/v0.9.4.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.4.rst
@@ -1,7 +1,7 @@
 .. _whatsnew_0940:
 
-v0.9.4 (TBD)
-------------------------
+v0.9.4 (anticipated December 2022)
+----------------------------------
 
 Deprecations
 ~~~~~~~~~~~~
@@ -10,6 +10,9 @@ Deprecations
 Enhancements
 ~~~~~~~~~~~~
 * Multiple code style issues fixed that were reported by LGTM analysis. (:issue:`1275`, :pull:`1559`)
+* Added a direct IAM model :py:func:`pvlib.iam.schlick` which can be used with
+  :py:func:`~pvlib.iam.marion_diffuse`, and a diffuse IAM model
+  :py:func:`pvlib.iam.schlick_diffuse` (:pull:`1562`, :issue:`1564`)
 * Added a function to calculate one of GHI, DHI, and DNI from values of the other two.
   :py:func:`~pvlib.irradiance.complete_irradiance`
   (:issue:`1565`, :pull:`1567`)
@@ -46,5 +49,9 @@ Contributors
 * Christian Orner (:ghuser:`chrisorner`)
 * Saurabh Aneja (:ghuser:`spaneja`)
 * Marcus Boumans (:ghuser:`bowie2211`)
+* Yu Xie (:ghuser:`xieyupku`)
+* Anton Driesse (:ghuser:`adriesse`)
+* Cliff Hansen (:ghuser:`cwhanse`)
+* Kevin Anderson (:ghuser:`kanderso-nrel`)
 * Karel De Brabandere (:ghuser:`kdebrab`)
 * Naman Priyadarshi (:ghuser:`Naman-Priyadarshi`)
diff --git a/pvlib/iam.py b/pvlib/iam.py
index a8592f4036..6e60b6f97d 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -541,7 +541,7 @@ def marion_diffuse(model, surface_tilt, **kwargs):
     ----------
     model : str
         The IAM function to evaluate across solid angle. Must be one of
-        `'ashrae', 'physical', 'martin_ruiz', 'sapm'`.
+        `'ashrae', 'physical', 'martin_ruiz', 'sapm', 'schlick'`.
 
     surface_tilt : numeric
         Surface tilt angles in decimal degrees.
@@ -592,6 +592,7 @@ def marion_diffuse(model, surface_tilt, **kwargs):
         'ashrae': ashrae,
         'sapm': sapm,
         'martin_ruiz': martin_ruiz,
+        'schlick': schlick,
     }
 
     try:
@@ -748,3 +749,123 @@ def marion_integrate(function, surface_tilt, region, num=None):
         Fd = pd.Series(Fd, surface_tilt.index)
 
     return Fd
+
+
+def schlick(aoi):
+    """
+    Determine incidence angle modifier (IAM) for direct irradiance using the
+    Schlick approximation to the Fresnel equations.
+
+    The Schlick approximation was proposed in [1]_ as a computationally
+    efficient alternative to computing the Fresnel factor in computer
+    graphics contexts.  This implementation is a normalized form of the
+    equation in [1]_ so that it can be used as a PV IAM model.
+    Unlike other IAM models, this model has no ability to describe
+    different reflection profiles.
+
+    In PV contexts, the Schlick approximation has been used as an analytically
+    integrable alternative to the Fresnel equations for estimating IAM
+    for diffuse irradiance [2]_.
+
+    Parameters
+    ----------
+    aoi : numeric
+        The angle of incidence (AOI) between the module normal vector and the
+        sun-beam vector. Angles of nan will result in nan. [degrees]
+
+    Returns
+    -------
+    iam : numeric
+        The incident angle modifier.
+
+    References
+    ----------
+    .. [1] Schlick, C. An inexpensive BRDF model for physically-based
+       rendering. Computer graphics forum 13 (1994).
+
+    .. [2] Xie, Y., M. Sengupta, A. Habte, A. Andreas, "The 'Fresnel Equations'
+       for Diffuse radiation on Inclined photovoltaic Surfaces (FEDIS)",
+       Renewable and Sustainable Energy Reviews, vol. 161, 112362. June 2022.
+       :doi:`10.1016/j.rser.2022.112362`
+
+    See Also
+    --------
+    pvlib.iam.schlick_diffuse
+    """
+    iam = 1 - (1 - cosd(aoi)) ** 5
+    iam = np.where(np.abs(aoi) >= 90.0, 0.0, iam)
+
+    # preserve input type
+    if np.isscalar(aoi):
+        iam = iam.item()
+    elif isinstance(aoi, pd.Series):
+        iam = pd.Series(iam, aoi.index)
+
+    return iam
+
+
+def schlick_diffuse(surface_tilt):
+    """
+    Determine the incidence angle modifiers (IAM) for diffuse sky and
+    ground-reflected irradiance on a tilted surface using the Schlick
+    incident angle model.
+
+    The diffuse iam values are calculated using an analytical integration
+    of the Schlick equation [1]_ over the portion of an isotropic sky and
+    isotropic foreground that is visible from the tilted surface [2]_.
+
+    Parameters
+    ----------
+    surface_tilt : numeric
+        Surface tilt angle measured from horizontal (e.g. surface facing
+        up = 0, surface facing horizon = 90). [degrees]
+
+    Returns
+    -------
+    iam_sky : numeric
+        The incident angle modifier for sky diffuse.
+
+    iam_ground : numeric
+        The incident angle modifier for ground-reflected diffuse.
+
+    References
+    ----------
+    .. [1] Schlick, C. An inexpensive BRDF model for physically-based
+       rendering. Computer graphics forum 13 (1994).
+
+    .. [2] Xie, Y., M. Sengupta, A. Habte, A. Andreas, "The 'Fresnel Equations'
+       for Diffuse radiation on Inclined photovoltaic Surfaces (FEDIS)",
+       Renewable and Sustainable Energy Reviews, vol. 161, 112362. June 2022.
+       :doi:`10.1016/j.rser.2022.112362`
+
+    See Also
+    --------
+    pvlib.iam.schlick
+    """
+    # these calculations are as in [2]_, but with the refractive index
+    # weighting coefficient w set to 1.0 (so it is omitted)
+
+    # relative transmittance of sky diffuse radiation by PV cover:
+    cosB = cosd(surface_tilt)
+    sinB = sind(surface_tilt)
+    cuk = (2 / (np.pi * (1 + cosB))) * (
+        (30/7)*np.pi - (160/21)*np.radians(surface_tilt) - (10/3)*np.pi*cosB
+        + (160/21)*cosB*sinB - (5/3)*np.pi*cosB*sinB**2 + (20/7)*cosB*sinB**3
+        - (5/16)*np.pi*cosB*sinB**4 + (16/105)*cosB*sinB**5
+    )  # Eq 4 in [2]
+
+    # relative transmittance of ground-reflected radiation by PV cover:
+    with np.errstate(divide='ignore', invalid='ignore'):  # Eq 6 in [2]
+        cug = 40 / (21 * (1 - cosB)) - (1 + cosB) / (1 - cosB) * cuk
+
+    cug = np.where(surface_tilt < 1e-6, 0, cug)
+
+    # respect input types:
+    if np.isscalar(surface_tilt):
+        cuk = cuk.item()
+        cug = cug.item()
+    elif isinstance(surface_tilt, pd.Series):
+        cuk = pd.Series(cuk, surface_tilt.index)
+        cug = pd.Series(cug, surface_tilt.index)
+
+    return cuk, cug
diff --git a/pvlib/tests/test_iam.py b/pvlib/tests/test_iam.py
index 4310ee837a..df4d9ee877 100644
--- a/pvlib/tests/test_iam.py
+++ b/pvlib/tests/test_iam.py
@@ -322,3 +322,48 @@ def test_marion_integrate_invalid():
 
     with pytest.raises(ValueError):
         _iam.marion_integrate(_iam.ashrae, 0, 'bad', 180)
+
+
+def test_schlick():
+    idx = pd.date_range('2019-01-01', freq='h', periods=9)
+    aoi = pd.Series([-180, -135, -90, -45, 0, 45, 90, 135, 180], idx)
+    expected = pd.Series([0, 0, 0, 0.99784451, 1.0, 0.99784451, 0, 0, 0], idx)
+
+    # scalars
+    for aoi_scalar, expected_scalar in zip(aoi, expected):
+        actual = _iam.schlick(aoi_scalar)
+        assert_allclose(expected_scalar, actual)
+
+    # numpy arrays
+    actual = _iam.schlick(aoi.values)
+    assert_allclose(expected.values, actual)
+
+    # pandas Series
+    actual = _iam.schlick(aoi)
+    assert_series_equal(expected, actual)
+
+
+def test_schlick_diffuse():
+    surface_tilt = np.array([0, 20, 70, 90])
+    # expected values calculated with marion_integrate and schlick
+    expected_sky = np.array([0.95238092, 0.96249934, 0.96228167, 0.95238094])
+    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)
+
+    # 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)
+
+    # 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,
+                        rtol=1e-6)