diff --git a/CHANGES.rst b/CHANGES.rst
index 664e60d0..90f96964 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,7 +1,19 @@
-1.6.1 (unreleased)
+1.7.0 (unreleased)
 ==================
 
+Changes to API
+--------------
+
+ramp_fitting
+~~~~~~~~~~~~
+
+- Add ``average_dark_current`` to calculations of poisson variance. [#243]
 
+Bug Fixes
+---------
+
+Other
+-----
 
 1.6.0 (2024-02-15)
 ==================
diff --git a/docs/stcal/ramp_fitting/description.rst b/docs/stcal/ramp_fitting/description.rst
index 37fc625c..9ea6830d 100644
--- a/docs/stcal/ramp_fitting/description.rst
+++ b/docs/stcal/ramp_fitting/description.rst
@@ -195,14 +195,16 @@ time in seconds (from the keyword TGROUP).
 The variance of the slope in a segment due to Poisson noise is:
 
 .. math::
-   var^P_{s} = \frac{ slope_{est} }{  tgroup \times gain\ (ngroups_{s} -1)}  \,,
+   var^P_{s} = \frac{ slope_{est} + darkcurrent}{  tgroup \times gain\ (ngroups_{s} -1)}  \,,
 
 where :math:`gain` is the gain for the pixel (from the GAIN reference file),
 in e/DN. The :math:`slope_{est}` is an overall estimated slope of the pixel,
 calculated by taking the median of the first differences of the groups that are
 unaffected by saturation and cosmic rays, in all integrations. This is a more
 robust estimate of the slope than the segment-specific slope, which may be noisy
-for short segments.
+for short segments. The contributions from the dark current are added when present;
+the value can be provided by the user during the `jwst.dark_current.DarkCurrentStep`,
+or it can be specified in scalar or 2D array form by the dark reference file.
 
 The combined variance of the slope of a segment is the sum of the variances:
 
diff --git a/src/stcal/ramp_fitting/ols_fit.py b/src/stcal/ramp_fitting/ols_fit.py
index f2c6d5a2..a52e2fba 100644
--- a/src/stcal/ramp_fitting/ols_fit.py
+++ b/src/stcal/ramp_fitting/ols_fit.py
@@ -542,8 +542,9 @@ def slice_ramp_data(ramp_data, start_row, nrows):
     err = ramp_data.err[:, :, start_row : start_row + nrows, :].copy()
     groupdq = ramp_data.groupdq[:, :, start_row : start_row + nrows, :].copy()
     pixeldq = ramp_data.pixeldq[start_row : start_row + nrows, :].copy()
+    average_dark_current = ramp_data.average_dark_current[start_row : start_row + nrows, :].copy()
 
-    ramp_data_slice.set_arrays(data, err, groupdq, pixeldq)
+    ramp_data_slice.set_arrays(data, err, groupdq, pixeldq, average_dark_current)
 
     if ramp_data.zeroframe is not None:
         ramp_data_slice.zeroframe = ramp_data.zeroframe[:, start_row : start_row + nrows, :].copy()
@@ -1138,7 +1139,8 @@ def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans)
             # Suppress harmless arithmetic warnings for now
             warnings.filterwarnings("ignore", ".*invalid value.*", RuntimeWarning)
             warnings.filterwarnings("ignore", ".*divide by zero.*", RuntimeWarning)
-            var_p4[num_int, :, rlo:rhi, :] = den_p3 * med_rates[rlo:rhi, :]
+            var_p4[num_int, :, rlo:rhi, :] = (den_p3 * med_rates[rlo:rhi, :]) + \
+                                             ramp_data.average_dark_current[rlo:rhi, :]
 
             # Find the segment variance due to read noise and convert back to DN
             var_r4[num_int, :, rlo:rhi, :] = num_r3 * den_r3 / gain_sect**2
@@ -1175,10 +1177,11 @@ def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans)
         # Huge variances correspond to non-existing segments, so are reset to 0
         #  to nullify their contribution.
         var_p3[var_p3 > utils.LARGE_VARIANCE_THRESHOLD] = 0.0
-        var_p3[:, med_rates <= 0.0] = 0.0
+        med_rate_mask = med_rates <= 0.0
+        var_p3[:, med_rate_mask] = 0.0
         warnings.resetwarnings()
 
-        var_p4[num_int, :, med_rates <= 0.0] = 0.0
+        var_p4[num_int, :, med_rate_mask] = ramp_data.average_dark_current[med_rate_mask][..., np.newaxis]
         var_both4[num_int, :, :, :] = var_r4[num_int, :, :, :] + var_p4[num_int, :, :, :]
         inv_var_both4[num_int, :, :, :] = 1.0 / var_both4[num_int, :, :, :]
 
diff --git a/src/stcal/ramp_fitting/ramp_fit.py b/src/stcal/ramp_fitting/ramp_fit.py
index 514429f1..11ba5b9e 100755
--- a/src/stcal/ramp_fitting/ramp_fit.py
+++ b/src/stcal/ramp_fitting/ramp_fit.py
@@ -53,10 +53,16 @@ def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False):
     """
     ramp_data = ramp_fit_class.RampData()
 
+    if not hasattr(model, 'average_dark_current'):
+        dark_current_array = np.zeros_like(model.pixeldq)
+    else:
+        dark_current_array = model.average_dark_current
+
     if isinstance(model.data, u.Quantity):
-        ramp_data.set_arrays(model.data.value, model.err.value, model.groupdq, model.pixeldq)
+        ramp_data.set_arrays(model.data.value, model.err.value, model.groupdq,
+                             model.pixeldq, dark_current_array)
     else:
-        ramp_data.set_arrays(model.data, model.err, model.groupdq, model.pixeldq)
+        ramp_data.set_arrays(model.data, model.err, model.groupdq, model.pixeldq, dark_current_array)
 
     # Attribute may not be supported by all pipelines.  Default is NoneType.
     drop_frames1 = model.meta.exposure.drop_frames1 if hasattr(model, "drop_frames1") else None
diff --git a/src/stcal/ramp_fitting/ramp_fit_class.py b/src/stcal/ramp_fitting/ramp_fit_class.py
index f8a78efd..ef2df6bc 100644
--- a/src/stcal/ramp_fitting/ramp_fit_class.py
+++ b/src/stcal/ramp_fitting/ramp_fit_class.py
@@ -6,6 +6,7 @@ def __init__(self):
         self.err = None
         self.groupdq = None
         self.pixeldq = None
+        self.average_dark_current = None
 
         # Meta information
         self.instrument_name = None
@@ -41,7 +42,7 @@ def __init__(self):
 
         self.current_integ = -1
 
-    def set_arrays(self, data, err, groupdq, pixeldq):
+    def set_arrays(self, data, err, groupdq, pixeldq, average_dark_current):
         """
         Set the arrays needed for ramp fitting.
 
@@ -62,12 +63,17 @@ def set_arrays(self, data, err, groupdq, pixeldq):
         pixeldq : ndarray (uint32)
             4-D array containing the pixel data quality information.  It has dimensions
             (nintegrations, ngroups, nrows, ncols)
+
+        average_dark_current : ndarray (float32)
+            2-D array containing the average dark current. It has
+            dimensions (nrows, ncols)
         """
         # Get arrays from the data model
         self.data = data
         self.err = err
         self.groupdq = groupdq
         self.pixeldq = pixeldq
+        self.average_dark_current = average_dark_current
 
     def set_meta(self, name, frame_time, group_time, groupgap, nframes, drop_frames1=None):
         """
diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py
index d8e90610..aab91a12 100644
--- a/tests/test_ramp_fitting.py
+++ b/tests/test_ramp_fitting.py
@@ -322,12 +322,13 @@ def jp_2326_test_setup():
     gdq = np.zeros((nints, ngroups, nrows, ncols), dtype=np.uint8)
     err = np.zeros((nints, ngroups, nrows, ncols))
     pdq = np.zeros((nrows, ncols), dtype=np.uint32)
+    dark_current = np.zeros((nrows, ncols))
 
     data[0, :, 0, 0] = ramp.copy()
     gdq[0, :, 0, 0] = dq.copy()
 
     ramp_data = RampData()
-    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pdq)
+    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pdq, average_dark_current=dark_current)
     ramp_data.set_meta(
         name="MIRI", frame_time=2.77504, group_time=2.77504, groupgap=0, nframes=1, drop_frames1=None
     )
@@ -417,6 +418,7 @@ def test_2_group_cases():
     rnoise = np.ones((1, npix)) * rnoise_val
     gain = np.ones((1, npix)) * gain_val
     pixeldq = np.zeros((1, npix), dtype=np.uint32)
+    dark_current = np.zeros((nrows, ncols), dtype=np.float32)
 
     data = np.zeros(dims, dtype=np.float32)  # Science data
     for k in range(npix):
@@ -433,7 +435,7 @@ def test_2_group_cases():
     # Setup the RampData class to run ramp fitting on.
     ramp_data = RampData()
 
-    ramp_data.set_arrays(data, err, groupdq, pixeldq)
+    ramp_data.set_arrays(data, err, groupdq, pixeldq, average_dark_current=dark_current)
 
     ramp_data.set_meta(
         name="NIRSPEC", frame_time=14.58889, group_time=14.58889, groupgap=0, nframes=1, drop_frames1=None
@@ -730,6 +732,8 @@ def create_zero_frame_data():
     pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32)
     gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8)
     zframe = np.ones(shape=(nints, nrows, ncols), dtype=np.float32)
+    dark_current = np.zeros((nrows, ncols), dtype=np.float32)
+
 
     # Create base ramps for each pixel in each integration.
     base_slope = 2000.0
@@ -756,7 +760,7 @@ def create_zero_frame_data():
 
     # Create RampData for testing.
     ramp_data = RampData()
-    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq)
+    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current)
     ramp_data.set_meta(
         name="NIRCam",
         frame_time=frame_time,
@@ -855,6 +859,8 @@ def create_only_good_0th_group_data():
     err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32)
     pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32)
     gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8)
+    dark_current = np.zeros((nrows, ncols), dtype=np.float32)
+
 
     # Create base ramps for each pixel in each integration.
     base_slope = 2000.0
@@ -877,7 +883,7 @@ def create_only_good_0th_group_data():
 
     # Create RampData for testing.
     ramp_data = RampData()
-    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq)
+    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current)
     ramp_data.set_meta(
         name="NIRCam",
         frame_time=frame_time,
@@ -1422,9 +1428,10 @@ def create_blank_ramp_data(dims, var, tm):
     err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32)
     pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32)
     gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8)
+    dark_current = np.zeros(shape=(nrows, ncols), dtype = np.float32)
 
     ramp_data = RampData()
-    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq)
+    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current)
     ramp_data.set_meta(
         name="NIRSpec",
         frame_time=frame_time,
@@ -1476,6 +1483,7 @@ def setup_inputs(dims, var, tm):
     err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32)
     pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32)
     gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8)
+    dark_current = np.zeros(shape=(nrows, ncols), dtype=np.float32)
 
     base_array = np.array([k + 1 for k in range(ngroups)])
     base, inc = 1.5, 1.5
@@ -1488,7 +1496,7 @@ def setup_inputs(dims, var, tm):
         data[c_int, :, :, :] = data[0, :, :, :].copy()
 
     ramp_data = RampData()
-    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq)
+    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current)
     ramp_data.set_meta(
         name="MIRI", frame_time=dtime, group_time=gtime, groupgap=0, nframes=nframes, drop_frames1=None
     )
diff --git a/tests/test_ramp_fitting_cases.py b/tests/test_ramp_fitting_cases.py
index 675e6a75..9e32813b 100644
--- a/tests/test_ramp_fitting_cases.py
+++ b/tests/test_ramp_fitting_cases.py
@@ -786,9 +786,10 @@ def create_blank_ramp_data(dims, var, timing, ts_name="NIRSpec"):
     err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32)
     pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32)
     gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8)
+    dark_current = np.zeros(shape=(nrows, ncols), dtype=np.float32)
 
     ramp_data = RampData()
-    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq)
+    ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current)
     ramp_data.set_meta(
         name=ts_name,
         frame_time=frame_time,
diff --git a/tests/test_ramp_fitting_gls_fit.py b/tests/test_ramp_fitting_gls_fit.py
index e940bd94..fef7c0a1 100644
--- a/tests/test_ramp_fitting_gls_fit.py
+++ b/tests/test_ramp_fitting_gls_fit.py
@@ -63,9 +63,11 @@ def setup_inputs(dims, gain, rnoise, group_time, frame_time):
     err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32)
     groupdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8)
     pixeldq = np.zeros(shape=(nrows, ncols), dtype=np.uint32)
+    dark_current = np.zeros(shape=(nrows, ncols), dtype=np.float32)
+
 
     # Set clas arrays
-    ramp_class.set_arrays(data, err, groupdq, pixeldq)
+    ramp_class.set_arrays(data, err, groupdq, pixeldq, average_dark_current=dark_current)
 
     # Set class meta
     ramp_class.set_meta(