diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 6e500b79..6e2a6574 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -5,7 +5,7 @@ build:
   apt_packages:
     - graphviz
   tools:
-    python: "3.13"
+    python: "3.14"
 
 sphinx:
   builder: html
diff --git a/CHANGES.rst b/CHANGES.rst
index 06eb977a..91103a6c 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -6,6 +6,14 @@ New Features
 
 - Added uncertainty propagation to ``specreduce.extract.BoxcarExtract`` and
   ``specreduce.extract.HorneExtract``. The extracted spectra have now proper uncertainties.
+- Added uncertainty propagation to ``specreduce.background.Background``. The
+  ``bkg_image()``, ``bkg_spectrum()``, ``sub_image()``, and ``sub_spectrum()`` methods
+  now return spectra with proper uncertainties. When input image has uncertainty, it is
+  propagated using variance formulas appropriate for the chosen statistic. When no
+  uncertainty is provided, it is estimated from the flux values in the background region.
+- Added optional ``sigma`` parameter to ``specreduce.background.Background`` for sigma
+  clipping outlier rejection in background estimation. Default is 5.0; set to ``None``
+  to disable.
 
 Other changes
 ^^^^^^^^^^^^^
diff --git a/docs/background.rst b/docs/background.rst
index 52747d4a..ab7a667a 100644
--- a/docs/background.rst
+++ b/docs/background.rst
@@ -31,4 +31,41 @@ and the background-subtracted image via `~specreduce.background.Background.sub_i
 (or ``image - bg``).
 
 The background and trace steps can be done iteratively, to refine an automated
-trace using the background-subtracted image as input.
\ No newline at end of file
+trace using the background-subtracted image as input.
+
+Outlier rejection
+-----------------
+
+The background estimation supports sigma clipping for outlier rejection, which is useful
+for removing cosmic rays and other artifacts from the background region. The ``sigma``
+parameter controls the number of standard deviations used for clipping (default is 5.0).
+Set ``sigma=None`` to disable sigma clipping.
+
+.. code-block:: python
+
+  # Use tighter sigma clipping for aggressive outlier rejection
+  bg = Background.two_sided(image, trace, separation=5, width=4, sigma=3.0)
+
+  # Disable sigma clipping
+  bg = Background.two_sided(image, trace, separation=5, width=4, sigma=None)
+
+Uncertainty propagation
+-----------------------
+
+The `~specreduce.background.Background` class propagates uncertainties through
+background subtraction. When the input image is an `~astropy.nddata.NDData` object
+with ``image.uncertainty`` provided, the uncertainties are propagated using variance
+formulas appropriate for the chosen statistic (``average`` or ``median``). When no
+uncertainty is provided, it is estimated from the variance of flux values in the
+background region.
+
+The background image, background spectrum, and background-subtracted outputs all
+include proper uncertainties:
+
+.. code-block:: python
+
+  bg = Background.two_sided(image, trace, separation=5, width=4)
+
+  # Access background with uncertainty
+  bkg_spec = bg.bkg_spectrum()
+  print(bkg_spec.uncertainty)
diff --git a/docs/conf.py b/docs/conf.py
index 1dc57bb1..970944d3 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -207,6 +207,17 @@
 # ArrayLike from numpy.typing is not properly resolved by intersphinx
 nitpick_ignore = [
     ("py:class", "ArrayLike"),
+    ("py:class", "numpy._typing.ArrayLike"),
+    ("py:class", "specutils.spectra.spectrum1d.Spectrum1D"),
+]
+
+# Ignore complex type annotations that can't be cross-referenced
+# These come from autodoc processing dataclass fields with Literal, Union, and generic types
+nitpick_ignore_regex = [
+    (r"py:obj", r"typing\.Literal\[.*"),      # Literal types with string values
+    (r"py:class", r"tuple\[.*"),               # Generic tuple types
+    (r"py:class", r"None \| dict\[.*"),        # Union types with dict
+    (r"py:class", r"dict\[.*"),                # Generic dict types
 ]
 
 # -- Options for linkcheck output -------------------------------------------
diff --git a/docs/extraction.rst b/docs/extraction.rst
index 789d12b8..4bd15b8e 100644
--- a/docs/extraction.rst
+++ b/docs/extraction.rst
@@ -44,6 +44,33 @@ flux is scaled according to the effective number of unmasked pixels. When using
 options (``filter`` or ``omit``), the non-finite values may be propagated or treated
 differently as documented in the API.
 
+Uncertainty propagation
+-----------------------
+
+Both `~specreduce.extract.BoxcarExtract` and `~specreduce.extract.HorneExtract`
+propagate uncertainties from the input image to the extracted 1D spectrum.
+
+For the input image, uncertainties can be provided in two ways:
+
+1. As part of an `~astropy.nddata.NDData` object via ``image.uncertainty``
+2. Via the ``variance`` parameter (for HorneExtract)
+
+The extracted spectrum includes the propagated uncertainty, which can be accessed via
+the ``uncertainty`` attribute:
+
+.. code-block:: python
+
+  extract = BoxcarExtract(image, trace, width=5)
+  spectrum = extract.spectrum
+  print(spectrum.uncertainty)
+
+For `~specreduce.extract.BoxcarExtract`, the uncertainty is propagated through the
+weighted sum over the extraction aperture. For `~specreduce.extract.HorneExtract`,
+the optimal extraction algorithm naturally produces properly weighted uncertainties.
+
+Calling the extraction methods
+------------------------------
+
 The previous examples in this section show how to initialize the BoxcarExtract
 or HorneExtract objects with their required parameters. To extract the 1D
 spectrum
diff --git a/environment-dev.yml b/environment-dev.yml
index cd4f65f5..7ffb956e 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -12,4 +12,8 @@ dependencies:
   - photutils
   - tox
   - sphinx-astropy
+  - sphinx-copybutton
+  - sphinx-design
+  - nbsphinx
+  - pydata-sphinx-theme
   - synphot
diff --git a/specreduce/background.py b/specreduce/background.py
index 2288fb31..2630cd49 100644
--- a/specreduce/background.py
+++ b/specreduce/background.py
@@ -5,6 +5,8 @@
 
 import numpy as np
 from astropy import units as u
+from astropy.nddata import VarianceUncertainty
+from astropy.stats import sigma_clip
 from astropy.utils.decorators import deprecated_attribute
 
 from specreduce.compat import SPECUTILS_LT_2, Spectrum
@@ -30,7 +32,7 @@ class Background(_ImageParser):
     ----------
     image
         Image with 2-D spectral image data
-    traces : List, `specreduce.tracing.Trace`, int, float
+    traces
         Individual or list of trace object(s) (or integers/floats to define
         FlatTraces) to extract the background. If None, a ``FlatTrace`` at the
         center of the image (according to ``disp_axis``) will be used.
@@ -57,9 +59,13 @@ class Background(_ImageParser):
             and the mask is dropped.
         - ``nan_fill``:  Pixels that are either masked or non-finite are replaced with nan,\
             and the mask is dropped.
-        - ``apply_mask_only``: The  image and mask are left unmodified.
-        - ``apply_nan_only``: The  image is left unmodified, the old mask is dropped, and a\
+        - ``apply_mask_only``: The image and mask are left unmodified.
+        - ``apply_nan_only``: The image is left unmodified, the old mask is dropped, and a\
             new mask is created based on non-finite values.
+    sigma
+        Number of standard deviations for sigma clipping outlier rejection.
+        Clipping is applied column-by-column along the cross-dispersion axis
+        within the background aperture. If None, no sigma clipping is performed.
 
 """
 
@@ -74,6 +80,7 @@ class Background(_ImageParser):
     disp_axis: int = 1
     crossdisp_axis: int = 0
     mask_treatment: MaskingOption = "apply"
+    sigma: None | float = 5.0
     _valid_mask_treatment_methods = (
         "apply",
         "ignore",
@@ -101,6 +108,13 @@ def __post_init__(self):
             raise ValueError("width must be positive")
         if self.width == 0:
             self._bkg_array = np.zeros(self.image.shape[self.disp_axis])
+            self._bkg_variance = np.zeros(self.image.shape[self.disp_axis])
+            self._variance_unit = self.image.unit**2
+            self._orig_uncty_type = (
+                type(self.image.uncertainty) if self.image.uncertainty is not None
+                else VarianceUncertainty
+            )
+            self._outlier_mask = np.zeros(self.image.shape, dtype=bool)
             return
 
         self._set_traces()
@@ -146,17 +160,39 @@ def __post_init__(self):
 
         self.bkg_wimage = bkg_wimage
 
+        # Apply sigma clipping for outlier rejection if enabled
+        if self.sigma is not None:
+            # Temporarily set mask to exclude pixels outside the background window for clipping
+            original_mask = img.mask.copy()
+            temp_mask = (self.bkg_wimage == 0) | original_mask
+            img.mask = temp_mask
+            # Apply sigma clipping column by column
+            clipped = sigma_clip(img, sigma=self.sigma, axis=self.crossdisp_axis,
+                                 masked=True, copy=False)
+            # Store only the NEW clipped mask (exclude pixels that were already masked)
+            self._outlier_mask = clipped.mask & ~temp_mask
+            # Restore original mask
+            img.mask = original_mask
+        else:
+            self._outlier_mask = np.zeros_like(img.data, dtype=bool)
+
+        combined_mask = img.mask | self._outlier_mask
+
         if self.statistic == "average":
-            self._bkg_array = np.ma.average(img, weights=self.bkg_wimage, axis=self.crossdisp_axis)
+            img_masked = np.ma.masked_array(img.data, combined_mask)
+            self._bkg_array = np.ma.average(img_masked, weights=self.bkg_wimage,
+                                            axis=self.crossdisp_axis)
 
         elif self.statistic == "median":
-            # combine where background weight image is 0 with image masked (which already
-            # accounts for non-finite data that wasn't already masked)
-            img.mask = np.logical_or(self.bkg_wimage == 0, self.image.mask)
-            self._bkg_array = np.ma.median(img, axis=self.crossdisp_axis)
+            combined_mask = combined_mask | (self.bkg_wimage == 0)
+            img_masked = np.ma.masked_array(img.data, combined_mask)
+            self._bkg_array = np.ma.median(img_masked, axis=self.crossdisp_axis)
         else:
             raise ValueError("statistic must be 'average' or 'median'")
 
+        # Compute background variance for uncertainty propagation
+        self._compute_bkg_variance(img)
+
     def _set_traces(self):
         """Determine `traces` from input. If an integer/float or list if int/float
         is passed in, use these to construct FlatTrace objects. These values
@@ -167,7 +203,7 @@ def _set_traces(self):
 
         if self.traces == []:
             # assume a flat trace at the image center if nothing is passed in.
-            trace_pos = self.image.shape[self.disp_axis] / 2.0
+            trace_pos = self.image.shape[self.crossdisp_axis] / 2.0
             self.traces = [FlatTrace(self.image, trace_pos)]
 
         if isinstance(self.traces, Trace):
@@ -192,13 +228,64 @@ def _set_traces(self):
                     "the middle of the image."
                 )
 
+    def _compute_bkg_variance(self, img):
+        """Compute background variance for uncertainty propagation.
+
+        Parameters
+        ----------
+        img : np.ma.MaskedArray
+            The masked image array used for background computation.
+        """
+        # Combined mask includes original mask and sigma-clipped pixels
+        combined_mask = img.mask | self._outlier_mask
+
+        # Get input variance (estimate from flux if not provided)
+        self._variance_unit = self.image.unit**2
+        if self.image.uncertainty is not None:
+            self._orig_uncty_type = type(self.image.uncertainty)
+            if self._orig_uncty_type == VarianceUncertainty:
+                var = self.image.uncertainty.array
+            else:
+                var = self.image.uncertainty.represent_as(VarianceUncertainty).array
+        else:
+            # Estimate variance from flux values in background region
+            self._orig_uncty_type = VarianceUncertainty
+            bkg_mask = self.bkg_wimage > 0
+            valid_mask = bkg_mask & ~combined_mask
+
+            # Compute sample variance along cross-dispersion axis for each column
+            masked_flux = np.ma.masked_array(self.image.data, ~valid_mask)
+            sample_variance = np.ma.var(masked_flux, axis=self.crossdisp_axis)
+
+            # Tile to full image shape (same variance for all pixels in a column)
+            var = np.tile(sample_variance.data, (self.image.shape[0], 1))
+
+        # Create masked variance array using combined mask
+        masked_variance = np.ma.masked_array(var, combined_mask)
+
+        if self.statistic == "average":
+            # Var(weighted_mean) = sum(w^2 * var) / (sum w)^2
+            weights_squared = self.bkg_wimage**2
+            numerator = np.ma.sum(weights_squared * masked_variance, axis=self.crossdisp_axis)
+            denominator = np.ma.sum(self.bkg_wimage, axis=self.crossdisp_axis) ** 2
+            self._bkg_variance = (numerator / denominator).data
+
+        elif self.statistic == "median":
+            # Var(median) ~ (pi/2) * var / n for approximately normal data
+            # Use mean variance of included pixels
+            bkg_mask = self.bkg_wimage > 0
+            valid_mask = bkg_mask & ~combined_mask
+            n_pixels = np.sum(valid_mask, axis=self.crossdisp_axis)
+            variance_sum = np.sum(np.where(valid_mask, var, 0), axis=self.crossdisp_axis)
+            mean_variance = variance_sum / n_pixels
+            self._bkg_variance = (np.pi / 2) * mean_variance / n_pixels
+
     @classmethod
     def two_sided(cls, image, trace_object, separation, **kwargs):
         """
         Determine the background from an image for subtraction centered around
         an input trace.
 
-
         Example: ::
 
             trace = FitTrace(image, guess=trace_pos)
@@ -206,37 +293,24 @@ def two_sided(cls, image, trace_object, separation, **kwargs):
 
         Parameters
         ----------
-        image : `~astropy.nddata.NDData`-like or array-like
+        image
             Image with 2-D spectral image data. Assumes cross-dispersion
             (spatial) direction is axis 0 and dispersion (wavelength)
             direction is axis 1.
-        trace_object : `~specreduce.tracing.Trace`
-            estimated trace of the spectrum to center the background traces
-        separation : float
-            separation from ``trace_object`` for the background regions
-        width : float
-            width of each background aperture in pixels
-        statistic : string
-            statistic to use when computing the background.  'average' will
-            account for partial pixel weights, 'median' will include all partial
-            pixels.
-        disp_axis : int
-            dispersion axis
-        crossdisp_axis : int
-            cross-dispersion axis
-        mask_treatment : string
-            The method for handling masked or non-finite data. Choice of ``filter``,
-            ``omit`, or ``zero_fill``. If `filter` is chosen, masked/non-finite data
-            will be filtered during the fit to each bin/column (along disp. axis) to
-            find the peak. If ``omit`` is chosen, columns along disp_axis with any
-            masked/non-finite data values will be fully masked (i.e, 2D mask is
-            collapsed to 1D and applied). If ``zero_fill`` is chosen, masked/non-finite
-            data will be replaced with 0.0 in the input image, and the mask will then
-            be dropped. For all three options, the input mask (optional on input
-            NDData object) will be combined with a mask generated from any non-finite
-            values in the image data.
-        """
+        trace_object
+            Estimated trace of the spectrum to center the background traces.
+        separation
+            Separation from ``trace_object`` for the background regions.
+        **kwargs
+            Additional keyword arguments passed to the `Background` constructor.
+            See `Background` for available options including ``width``, ``statistic``,
+            ``disp_axis``, ``crossdisp_axis``, ``mask_treatment``, and ``sigma``.
 
+        Returns
+        -------
+        `Background`
+            A Background object with two-sided background regions.
+        """
         image = _ImageParser._get_data_from_image(image) if image is not None else cls.image
         kwargs["traces"] = [trace_object - separation, trace_object + separation]
         return cls(image=image, **kwargs)
@@ -254,42 +328,30 @@ def one_sided(cls, image, trace_object, separation, **kwargs):
 
         Parameters
         ----------
-        image : `~astropy.nddata.NDData`-like or array-like
+        image
             Image with 2-D spectral image data. Assumes cross-dispersion
             (spatial) direction is axis 0 and dispersion (wavelength)
             direction is axis 1.
-        trace_object : `~specreduce.tracing.Trace`
-            Estimated trace of the spectrum to center the background traces
-        separation : float
-            Separation from ``trace_object`` for the background, positive will be
+        trace_object
+            Estimated trace of the spectrum to center the background trace.
+        separation
+            Separation from ``trace_object`` for the background. Positive will be
             above the trace, negative below.
-        width : float
-            Width of each background aperture in pixels
-        statistic : string
-            Statistic to use when computing the background.  'average' will
-            account for partial pixel weights, 'median' will include all partial
-            pixels.
-        disp_axis : int
-            Dispersion axis
-        crossdisp_axis : int
-            Cross-dispersion axis
-        mask_treatment : string
-            The method for handling masked or non-finite data. Choice of ``filter``,
-            ``omit``, or ``zero_fill``. If `filter` is chosen, masked/non-finite data
-            will be filtered during the fit to each bin/column (along disp. axis) to
-            find the peak. If ``omit`` is chosen, columns along disp_axis with any
-            masked/non-finite data values will be fully masked (i.e, 2D mask is
-            collapsed to 1D and applied). If ``zero_fill`` is chosen, masked/non-finite
-            data will be replaced with 0.0 in the input image, and the mask will then
-            be dropped. For all three options, the input mask (optional on input
-            NDData object) will be combined with a mask generated from any non-finite
-            values in the image data.
+        **kwargs
+            Additional keyword arguments passed to the `Background` constructor.
+            See `Background` for available options including ``width``, ``statistic``,
+            ``disp_axis``, ``crossdisp_axis``, ``mask_treatment``, and ``sigma``.
+
+        Returns
+        -------
+        `Background`
+            A Background object with a one-sided background region.
         """
         image = _ImageParser._get_data_from_image(image) if image is not None else cls.image
         kwargs["traces"] = [trace_object + separation]
         return cls(image=image, **kwargs)
 
-    def bkg_image(self, image=None):
+    def bkg_image(self, image=None) -> Spectrum:
         """
         Expose the background tiled to the dimension of ``image``.
 
@@ -303,59 +365,75 @@ def bkg_image(self, image=None):
 
         Returns
         -------
-        spec : `~specutils.Spectrum1D`
-            Spectrum object with same shape as ``image``.
+        spec : `~specutils.Spectrum`
+            Spectrum object with same shape as ``image``, including uncertainty.
         """
         image = self._parse_image(image)
         arr = np.tile(self._bkg_array, (image.shape[0], 1))
+        var_arr = np.tile(self._bkg_variance, (image.shape[0], 1))
+        uncertainty = VarianceUncertainty(var_arr * self._variance_unit).represent_as(
+            self._orig_uncty_type
+        )
+
         if SPECUTILS_LT_2:
             kwargs = {}
         else:
             kwargs = {"spectral_axis_index": arr.ndim - 1}
         return Spectrum(
-            arr * image.unit,
-            spectral_axis=image.spectral_axis, **kwargs
+            arr * image.unit, spectral_axis=image.spectral_axis, uncertainty=uncertainty, **kwargs
         )
 
-    def bkg_spectrum(self, image=None, bkg_statistic=None):
+    def bkg_spectrum(self, image=None, bkg_statistic=None) -> Spectrum:
         """
         Expose the 1D spectrum of the background.
 
         Parameters
         ----------
-        image : `~astropy.nddata.NDData`-like or array-like, optional
-            Image with 2-D spectral image data. Assumes cross-dispersion
+        image
+            Image with 2D spectral image data. Assumes cross-dispersion
             (spatial) direction is axis 0 and dispersion (wavelength)
             direction is axis 1. If None, will extract the background
-            from ``image`` used to initialize the class. [default: None]
+            from ``image`` used to initialize the class.
+        bkg_statistic
+            Deprecated. Use ``statistic`` parameter in the `Background`
+            constructor instead.
 
         Returns
         -------
-        spec : `~specutils.Spectrum1D`
-            The background 1-D spectrum, with flux expressed in the same
-            units as the input image (or DN if none were provided).
+        `~specutils.Spectrum1D`
+            The background 1D spectrum, with flux and uncertainty expressed
+            in the same units as the input image (or DN if none were provided).
         """
         if bkg_statistic is not None:
-            warnings.warn("'bkg_statistic' is deprecated and will be removed in a future release. "
-                          "Please use the 'statistic' argument in the Background initializer instead.",  # noqa
-                          DeprecationWarning,)
+            warnings.warn(
+                "'bkg_statistic' is deprecated and will be removed in a future release. "
+                "Please use the 'statistic' argument in the Background initializer instead.",  # noqa
+                DeprecationWarning,
+            )
 
-        return Spectrum(self._bkg_array * self.image.unit, self.image.spectral_axis)
+        uncertainty = VarianceUncertainty(self._bkg_variance * self._variance_unit).represent_as(
+            self._orig_uncty_type
+        )
+        return Spectrum(
+            self._bkg_array * self.image.unit, self.image.spectral_axis, uncertainty=uncertainty
+        )
 
-    def sub_image(self, image=None):
+    def sub_image(self, image=None) -> Spectrum:
         """
         Subtract the computed background from image.
 
         Parameters
         ----------
-        image : nddata-compatible image or None
-            image with 2D spectral image data.  If None, will extract
-            the background from ``image`` used to initialize the class.
+        image
+            `~astropy.nddata.NDData`-like or array-like 2D spectral image data.
+            If None, will extract the background from ``image`` used to
+            initialize the class.
 
         Returns
         -------
-        spec : `~specutils.Spectrum`
-            Spectrum object with same shape as ``image``.
+        `~specutils.Spectrum`
+            Spectrum object with same shape as ``image``, including propagated
+            uncertainty.
         """
         image = self._parse_image(image)
 
@@ -370,7 +448,7 @@ def sub_image(self, image=None):
         # https://docs.astropy.org/en/stable/nddata/mixins/ndarithmetic.html
         return image.subtract(self.bkg_image(image), **kwargs)
 
-    def sub_spectrum(self, image=None):
+    def sub_spectrum(self, image=None) -> Spectrum:
         """
         Expose the 1D spectrum of the background-subtracted image.
 
@@ -383,22 +461,45 @@ def sub_spectrum(self, image=None):
         Returns
         -------
         spec : `~specutils.Spectrum`
-            The background 1D spectrum, with flux expressed in the same
-            units as the input image (or u.DN if none were provided) and
-            the spectral axis expressed in pixel units.
+            The background-subtracted 1D spectrum, with flux and uncertainty
+            expressed in the same units as the input image (or u.DN if none
+            were provided) and the spectral axis expressed in pixel units.
         """
-        sub_image = self.sub_image(image=image)
+        sub_img = self.sub_image(image=image)
 
         try:
-            return sub_image.collapse(np.nansum, axis=self.crossdisp_axis)
+            result = sub_img.collapse(np.nansum, axis=self.crossdisp_axis)
         except u.UnitTypeError:
             # can't collapse with a spectral axis in pixels because
             # SpectralCoord only allows frequency/wavelength equivalent units...
-            ext1d = np.nansum(sub_image.flux, axis=self.crossdisp_axis)
-            return Spectrum(ext1d, spectral_axis=sub_image.spectral_axis)
+            ext1d = np.nansum(sub_img.flux, axis=self.crossdisp_axis)
+            result = Spectrum(ext1d, spectral_axis=sub_img.spectral_axis)
+
+        # Propagate uncertainty: Var(sum) = sum of variances
+        # Must convert to variance before summing, then convert back to original type
+        if sub_img.uncertainty is not None:
+            var_uncert = sub_img.uncertainty.represent_as(VarianceUncertainty)
+            var_sum = np.nansum(var_uncert.array, axis=self.crossdisp_axis)
+            uncertainty = VarianceUncertainty(var_sum * var_uncert.unit).represent_as(
+                self._orig_uncty_type
+            )
+            result = Spectrum(
+                result.flux, spectral_axis=result.spectral_axis, uncertainty=uncertainty
+            )
+        return result
 
     def __rsub__(self, image):
         """
         Subtract the background from an image.
+
+        Parameters
+        ----------
+        image
+            `~astropy.nddata.NDData`-like or array-like 2D spectral image data.
+
+        Returns
+        -------
+        `~specutils.Spectrum`
+            Background-subtracted image with propagated uncertainty.
         """
         return self.sub_image(image)
diff --git a/specreduce/tests/test_background.py b/specreduce/tests/test_background.py
index 9dda5329..d4d99991 100644
--- a/specreduce/tests/test_background.py
+++ b/specreduce/tests/test_background.py
@@ -1,7 +1,7 @@
 import astropy.units as u
 import numpy as np
 import pytest
-from astropy.nddata import NDData
+from astropy.nddata import NDData, VarianceUncertainty, StdDevUncertainty, InverseVariance
 
 from specreduce.background import Background
 from specreduce.compat import Spectrum
@@ -131,7 +131,7 @@ def test_trace_inputs(mk_test_img_raw):
     # When `Background` object is created with no Trace object passed in it should
     # create a FlatTrace in the middle of the image (according to disp. axis)
     background = Background(image, width=5)
-    assert np.all(background.traces[0].trace.data == image.shape[1] / 2.0)
+    assert np.all(background.traces[0].trace.data == image.shape[0] / 2.0)
 
     # FlatTrace(s) should be created if number or list of numbers is passed in for `traces`
     background = Background(image, 10.0, width=5)
@@ -338,3 +338,364 @@ def test_sub_bkg_image(self):
 
         assert np.all(np.isfinite(subtracted_img_zero_fill.data))
         assert np.all(subtracted_img_zero_fill.mask == 0)
+
+
+@pytest.fixture
+def mk_bgk_img():
+    nrows, ncols = 10, 20
+    var = 4.0
+    img = Spectrum(
+        np.ones((nrows, ncols)) * u.DN,
+        uncertainty=VarianceUncertainty(np.full((nrows, ncols), var) * u.DN**2),
+        spectral_axis=np.arange(ncols) * u.pix
+    )
+    trace = FlatTrace(img, nrows // 2)
+    return trace, img, var, nrows, ncols
+
+
+def test_background_uncertainty_average(mk_bgk_img):
+    """Test background uncertainty estimation with 'average' statistic."""
+    trace, img, var, nrows, ncols = mk_bgk_img
+    bg = Background(img, trace, width=4, statistic="average")
+
+    # Check that _bkg_variance is computed
+    assert hasattr(bg, "_bkg_variance")
+    assert bg._bkg_variance is not None
+    assert len(bg._bkg_variance) == ncols
+
+    # Variance should be positive and less than input (averaging reduces variance)
+    assert np.all(bg._bkg_variance > 0)
+    assert np.all(bg._bkg_variance < var)
+
+    weights = bg.bkg_wimage
+    weights_sum = np.sum(weights, axis=0)
+    weights_sq_sum = np.sum(weights ** 2, axis=0)
+    expected_variance = (weights_sq_sum * var) / (weights_sum ** 2)
+    assert np.allclose(bg._bkg_variance, expected_variance)
+
+    # Check that bkg_spectrum has uncertainty
+    bkg_spec = bg.bkg_spectrum()
+    assert bkg_spec.uncertainty is not None
+    assert isinstance(bkg_spec.uncertainty, VarianceUncertainty)
+    assert np.allclose(bkg_spec.uncertainty.array, bg._bkg_variance)
+
+
+def test_background_uncertainty_median(mk_bgk_img):
+    """Test background uncertainty estimation with 'median' statistic."""
+    trace, img, var, nrows, ncols = mk_bgk_img
+    bg = Background(img, trace, width=4, statistic="median")
+
+    # Check that _bkg_variance is computed
+    assert hasattr(bg, "_bkg_variance")
+    assert bg._bkg_variance is not None
+
+    # Variance should be positive
+    assert np.all(bg._bkg_variance > 0)
+    assert np.all(np.isfinite(bg._bkg_variance))
+
+    n_pixels = np.sum(bg.bkg_wimage > 0, axis=0)
+    expected_variance = (np.pi / 2) * var / n_pixels
+    assert np.allclose(bg._bkg_variance, expected_variance, rtol=0.01)
+
+
+def test_bkg_image_has_uncertainty(mk_bgk_img):
+    """Test that bkg_image returns Spectrum with uncertainty."""
+    trace, img, var, nrows, ncols = mk_bgk_img
+    bg = Background(img, trace, width=4)
+    bkg_img = bg.bkg_image()
+
+    # Check uncertainty exists
+    assert bkg_img.uncertainty is not None
+    assert isinstance(bkg_img.uncertainty, VarianceUncertainty)
+
+    # Check shape matches image
+    assert bkg_img.uncertainty.array.shape == (nrows, ncols)
+
+    # Check values are tiled correctly (same value in each column)
+    for col in range(ncols):
+        assert np.allclose(bkg_img.uncertainty.array[:, col], bg._bkg_variance[col])
+
+
+def test_bkg_spectrum_has_uncertainty(mk_bgk_img):
+    """Test that bkg_spectrum returns Spectrum with uncertainty."""
+    trace, img, var, nrows, ncols = mk_bgk_img
+    bg = Background(img, trace, width=4)
+    bkg_spec = bg.bkg_spectrum()
+
+    # Check uncertainty exists and is 1D
+    assert bkg_spec.uncertainty is not None
+    assert isinstance(bkg_spec.uncertainty, VarianceUncertainty)
+    assert bkg_spec.uncertainty.array.ndim == 1
+    assert len(bkg_spec.uncertainty.array) == ncols
+
+
+def test_sub_image_propagates_uncertainty(mk_bgk_img):
+    """Test that sub_image propagates uncertainties correctly."""
+    trace, img, var, nrows, ncols = mk_bgk_img
+    bg = Background(img, trace, width=4)
+
+    sub_img = bg.sub_image()
+
+    # Check uncertainty exists
+    assert sub_img.uncertainty is not None
+    assert isinstance(sub_img.uncertainty, VarianceUncertainty)
+
+    # For subtraction: Var(A - B) = Var(A) + Var(B)
+    # Image variance is 4.0, background variance is computed
+    bkg_variance = bg._bkg_variance[0]  # Same for all columns
+    expected_variance = var + bkg_variance
+    assert np.allclose(sub_img.uncertainty.array, expected_variance, rtol=0.01)
+
+
+def test_sub_spectrum_propagates_uncertainty(mk_bgk_img):
+    """Test that sub_spectrum propagates uncertainties correctly."""
+    trace, img, var, nrows, ncols = mk_bgk_img
+    bg = Background(img, trace, width=4)
+
+    sub_spec = bg.sub_spectrum()
+
+    # Check uncertainty exists
+    assert sub_spec.uncertainty is not None
+    assert isinstance(sub_spec.uncertainty, VarianceUncertainty)
+
+    # Check it's 1D with correct length
+    assert sub_spec.uncertainty.array.ndim == 1
+    assert len(sub_spec.uncertainty.array) == ncols
+
+    # Variance should be positive and finite
+    assert np.all(np.isfinite(sub_spec.uncertainty.array))
+    assert np.all(sub_spec.uncertainty.array > 0)
+
+
+def test_background_uncertainty_with_mask():
+    """Test that masked pixels don't contribute to uncertainty calculation."""
+    nrows, ncols = 10, 20
+    flux = np.ones((nrows, ncols))
+    variance = np.full((nrows, ncols), 4.0)
+    mask = np.zeros((nrows, ncols), dtype=bool)
+
+    # Mask some pixels in the background region for first few columns
+    mask[3:7, 0:5] = True
+
+    img = Spectrum(
+        flux * u.DN,
+        uncertainty=VarianceUncertainty(variance * u.DN**2),
+        mask=mask,
+        spectral_axis=np.arange(ncols) * u.pix
+    )
+
+    trace = FlatTrace(img, nrows // 2)
+    bg = Background(img, trace, width=4)
+
+    # Uncertainty should be computed
+    assert bg._bkg_variance is not None
+    assert np.all(np.isfinite(bg._bkg_variance))
+
+    # Columns with masked pixels should have different (larger) variance
+    # because fewer pixels contribute to the average
+    # (This depends on whether masked pixels fall in the background window)
+
+
+def test_background_no_input_uncertainty():
+    """Test that background estimates variance from flux when no uncertainty provided."""
+    nrows, ncols = 10, 20
+
+    # Create image with known noise pattern (no uncertainty provided)
+    np.random.seed(42)
+    noise_stddev = 2.0
+    flux = 100.0 + np.random.normal(0, noise_stddev, (nrows, ncols))
+    img = Spectrum(
+        flux * u.DN,
+        spectral_axis=np.arange(ncols) * u.pix
+    )
+
+    trace = FlatTrace(img, nrows // 2)
+    bg = Background(img, trace, width=4)
+
+    # Should have variance estimated from flux
+    assert hasattr(bg, "_bkg_variance")
+    assert bg._bkg_variance is not None
+
+    # Variance should be positive and finite (estimated from flux, not defaulting to 1.0)
+    assert np.all(bg._bkg_variance >= 0)
+    assert np.all(np.isfinite(bg._bkg_variance))
+
+    # Verify variance was computed from flux by checking it varies with the data
+    # (if it defaulted to 1.0, all values would be identical)
+    # With noisy data, different columns should have different sample variances
+    assert not np.allclose(bg._bkg_variance, 1.0)  # Not defaulting to 1.0
+
+    # bkg_spectrum should have uncertainty
+    bkg_spec = bg.bkg_spectrum()
+    assert bkg_spec.uncertainty is not None
+
+
+def test_background_uncertainty_stddev_type():
+    """Test that StdDevUncertainty type is preserved in Background output."""
+    nrows, ncols = 10, 20
+    stddev = 2.0  # stddev=2 means variance=4
+    img = Spectrum(
+        np.ones((nrows, ncols)) * 10 * u.DN,
+        uncertainty=StdDevUncertainty(np.full((nrows, ncols), stddev) * u.DN),
+        spectral_axis=np.arange(ncols) * u.pix
+    )
+
+    trace = FlatTrace(img, nrows // 2)
+    bg = Background(img, trace, width=4)
+
+    # bkg_spectrum should return StdDevUncertainty
+    bkg_spec = bg.bkg_spectrum()
+    assert bkg_spec.uncertainty is not None
+    assert isinstance(bkg_spec.uncertainty, StdDevUncertainty)
+
+    # bkg_image should return StdDevUncertainty
+    bkg_img = bg.bkg_image()
+    assert bkg_img.uncertainty is not None
+    assert isinstance(bkg_img.uncertainty, StdDevUncertainty)
+
+    # sub_spectrum should return StdDevUncertainty
+    sub_spec = bg.sub_spectrum()
+    assert sub_spec.uncertainty is not None
+    assert isinstance(sub_spec.uncertainty, StdDevUncertainty)
+
+    # Values should be positive and finite
+    assert np.all(bkg_spec.uncertainty.array > 0)
+    assert np.all(np.isfinite(bkg_spec.uncertainty.array))
+
+
+def test_background_uncertainty_inverse_variance_type():
+    """Test that InverseVariance type is preserved in Background output."""
+    nrows, ncols = 10, 20
+    ivar = 0.25  # ivar=0.25 means variance=4
+    img = Spectrum(
+        np.ones((nrows, ncols)) * 10 * u.DN,
+        uncertainty=InverseVariance(np.full((nrows, ncols), ivar) / u.DN**2),
+        spectral_axis=np.arange(ncols) * u.pix
+    )
+
+    trace = FlatTrace(img, nrows // 2)
+    bg = Background(img, trace, width=4)
+
+    # bkg_spectrum should return InverseVariance
+    bkg_spec = bg.bkg_spectrum()
+    assert bkg_spec.uncertainty is not None
+    assert isinstance(bkg_spec.uncertainty, InverseVariance)
+
+    # bkg_image should return InverseVariance
+    bkg_img = bg.bkg_image()
+    assert bkg_img.uncertainty is not None
+    assert isinstance(bkg_img.uncertainty, InverseVariance)
+
+    # sub_spectrum should return InverseVariance
+    sub_spec = bg.sub_spectrum()
+    assert sub_spec.uncertainty is not None
+    assert isinstance(sub_spec.uncertainty, InverseVariance)
+
+    # Values should be positive and finite
+    assert np.all(bkg_spec.uncertainty.array > 0)
+    assert np.all(np.isfinite(bkg_spec.uncertainty.array))
+
+
+def test_sigma_clipping_rejects_outliers():
+    """Test that sigma clipping rejects outlier pixels in background region."""
+    nrows, ncols = 20, 30
+    true_bkg = 10.0
+
+    # Create image with uniform background
+    flux = np.ones((nrows, ncols)) * true_bkg
+
+    # Add extreme outliers in background region (should be clipped)
+    flux[2, 5] = 1000.0  # extreme outlier
+    flux[17, 10] = -500.0  # extreme outlier
+
+    img = Spectrum(
+        flux * u.DN,
+        uncertainty=VarianceUncertainty(np.ones((nrows, ncols)) * u.DN**2),
+        spectral_axis=np.arange(ncols) * u.pix
+    )
+
+    trace = FlatTrace(img, nrows // 2)
+
+    # With sigma clipping enabled (low sigma to ensure clipping)
+    bg_clipped = Background(img, trace, width=4, sigma=3.0)
+
+    # Without sigma clipping
+    bg_no_clip = Background(img, trace, width=4, sigma=None)
+
+    # Sigma clip mask should mark the outliers (if they fall in bkg region)
+    assert hasattr(bg_clipped, "_outlier_mask")
+    assert bg_clipped._outlier_mask.shape == img.data.shape
+
+    # Background with clipping should be closer to true value
+    # than without clipping (for columns with outliers)
+    bkg_clipped = bg_clipped.bkg_spectrum().flux.value
+    bkg_no_clip = bg_no_clip.bkg_spectrum().flux.value
+
+    # Column 5 has outlier at row 2 - clipped version should be closer to true_bkg
+    if bg_clipped._outlier_mask[2, 5]:  # outlier was in bkg region and clipped
+        assert np.abs(bkg_clipped[5] - true_bkg) < np.abs(bkg_no_clip[5] - true_bkg)
+
+
+def test_sigma_clipping_disabled():
+    """Test that sigma=None disables sigma clipping."""
+    nrows, ncols = 10, 20
+    flux = np.ones((nrows, ncols)) * 10.0
+    flux[2, 5] = 100.0
+
+    img = Spectrum(
+        flux * u.DN,
+        spectral_axis=np.arange(ncols) * u.pix
+    )
+
+    trace = FlatTrace(img, nrows // 2)
+    bg = Background(img, trace, width=4, sigma=None)
+
+    # Sigma clip mask should be all False
+    assert hasattr(bg, "_outlier_mask")
+    assert not np.any(bg._outlier_mask)
+
+
+def test_sigma_clipping_preserves_image_mask():
+    """Test that sigma clipping doesn't modify the original image mask."""
+    nrows, ncols = 10, 20
+    flux = np.ones((nrows, ncols)) * 10.0
+    flux[2, 5] = 1000.0  # extreme outlier
+
+    # Create image with existing mask
+    original_mask = np.zeros((nrows, ncols), dtype=bool)
+    original_mask[3, 7] = True
+    original_mask[5, 12] = True
+
+    img = Spectrum(
+        flux * u.DN,
+        mask=original_mask.copy(),
+        spectral_axis=np.arange(ncols) * u.pix
+    )
+
+    trace = FlatTrace(img, nrows // 2)
+    bg = Background(img, trace, width=4, sigma=3.0)
+
+    # Original image mask should be unchanged
+    assert np.array_equal(bg.image.mask, original_mask)
+
+    # Sigma clip mask should be stored separately
+    assert hasattr(bg, "_outlier_mask")
+    assert not np.array_equal(bg._outlier_mask, original_mask)
+
+
+def test_sigma_clipping_default():
+    """Test that default sigma=5 is used when not specified."""
+    nrows, ncols = 10, 20
+    flux = np.ones((nrows, ncols)) * 10.0
+
+    img = Spectrum(
+        flux * u.DN,
+        spectral_axis=np.arange(ncols) * u.pix
+    )
+
+    trace = FlatTrace(img, nrows // 2)
+    bg = Background(img, trace, width=4)
+
+    # Default sigma should be 5
+    assert bg.sigma == 5.0
+    assert hasattr(bg, "_outlier_mask")