diff --git a/specutils/io/default_loaders/sdss_v.py b/specutils/io/default_loaders/sdss_v.py
index d3575f3d9..3aaec2250 100644
--- a/specutils/io/default_loaders/sdss_v.py
+++ b/specutils/io/default_loaders/sdss_v.py
@@ -481,7 +481,7 @@ def load_sdss_spec_list(file_obj, **kwargs):
         return SpectrumList(spectra)
 
 
-def _load_BOSS_HDU(hdulist: HDUList, hdu: int, **kwargs):
+def _load_BOSS_HDU(hdulist: HDUList, hdu: int, model: bool = False, **kwargs):
     """
     HDU processor for BOSS spectra redux HDU's
 
@@ -503,7 +503,14 @@ def _load_BOSS_HDU(hdulist: HDUList, hdu: int, **kwargs):
     flux_unit = Unit("1e-17 erg / (Angstrom cm2 s)")  # NOTE: hardcoded unit
     spectral_axis = Quantity(10**hdulist[hdu].data["LOGLAM"], unit=Angstrom)
 
-    flux = Quantity(hdulist[hdu].data["FLUX"], unit=flux_unit)
+    if model and 'MODEL' not in hdulist[hdu].data.names:
+        raise ValueError(f"MODEL column not found in HDU{hdu}. Cannot load model spectrum.")
+
+    if model:
+        # set the model spectrum as the "flux" ; until specutils supports models
+        flux = Quantity(hdulist[hdu].data["MODEL"], unit=flux_unit)
+    else:
+        flux = Quantity(hdulist[hdu].data["FLUX"], unit=flux_unit)
     # no e_flux, so we use inverse of variance
     ivar = InverseVariance(hdulist[hdu].data["IVAR"])
 
@@ -522,6 +529,7 @@ def _load_BOSS_HDU(hdulist: HDUList, hdu: int, **kwargs):
     meta = dict()
     meta["header"] = hdulist[0].header
     meta["name"] = hdulist[hdu].name
+    meta["is_model"] = model
 
     return Spectrum(spectral_axis=spectral_axis,
                       flux=flux,
diff --git a/specutils/io/default_loaders/tests/test_sdss_v.py b/specutils/io/default_loaders/tests/test_sdss_v.py
index c14cf3568..8edf3b192 100644
--- a/specutils/io/default_loaders/tests/test_sdss_v.py
+++ b/specutils/io/default_loaders/tests/test_sdss_v.py
@@ -420,7 +420,7 @@ def apVisit_HDUList():
     return hdulist
 
 
-def spec_HDUList(n_spectra):
+def spec_HDUList(n_spectra, model=True):
     """Mock an BOSS spec HDUList of n_spectra spectra + 1 coadd."""
     np.random.seed(20)
 
@@ -437,10 +437,7 @@ def spec_HDUList(n_spectra):
     hdulist = fits.HDUList()
     hdulist.append(fits.PrimaryHDU(header=hdr))
 
-    # Init the key HDU's (flux, error, bitmask, spectral)
-    names = ["COADD", "SPALL", "ZALL", "ZLINE"]
-    for i in range(4):
-        hdu = fits.BinTableHDU.from_columns([
+    cols = [
             fits.Column(name="FLUX", format="E", array=np.random.random(10)),
             fits.Column(name="LOGLAM",
                         format="E",
@@ -451,22 +448,18 @@ def spec_HDUList(n_spectra):
                         array=np.random.random(10)),
             fits.Column(name="OR_MASK", format="E",
                         array=np.random.random(10)),
-        ])
+        ]
+    if model:
+        cols.append(fits.Column(name="MODEL", format="E", array=np.random.random(10)))
+
+    # Init the key HDU's (flux, error, bitmask, spectral)
+    names = ["COADD", "SPALL", "ZALL", "ZLINE"]
+    for i in range(4):
+        hdu = fits.BinTableHDU.from_columns(cols)
         hdu.name = names[i]
         hdulist.append(hdu)
     for i in range(n_spectra):
-        hdu = fits.BinTableHDU.from_columns([
-            fits.Column(name="LOGLAM",
-                        format="E",
-                        array=np.random.random(10).sort()),
-            fits.Column(name="FLUX", format="E", array=np.random.random(10)),
-            fits.Column(name="IVAR", format="E", array=np.random.random(10)),
-            fits.Column(name="AND_MASK",
-                        format="E",
-                        array=np.random.random(10)),
-            fits.Column(name="OR_MASK", format="E",
-                        array=np.random.random(10)),
-        ])
+        hdu = fits.BinTableHDU.from_columns(cols)
         hdu.name = f"spectrum{i}"
         hdulist.append(hdu)
 
@@ -937,6 +930,28 @@ def test_spec_1d_fail_hdu(file_obj, hdu):
     os.remove(tmpfile)
 
 
+def test_spec_model():
+    """Test if model spectrum can be loaded"""
+    tmpfile = "spec-temp.fits"
+    spec_HDUList(1, model=True).writeto(tmpfile, overwrite=True)
+
+    flux = Spectrum.read(tmpfile, hdu=1, model=False)
+    model = Spectrum.read(tmpfile, hdu=1, model=True)
+    assert flux != model
+    assert flux.meta["is_model"] is False
+    assert model.meta["is_model"] is True
+    os.remove(tmpfile)
+
+def test_spec_model_fail():
+    """Test if model spectrum is not available """
+    tmpfile = "spec-temp.fits"
+    spec_HDUList(1, model=False).writeto(tmpfile, overwrite=True)
+
+    with pytest.raises(ValueError, match='MODEL column not found in HDU'):
+        Spectrum.read(tmpfile, hdu=1, model=True)
+    os.remove(tmpfile)
+
+
 @pytest.mark.parametrize(
     "file_obj,idx",
     [