diff --git a/autotest/t003_test.py b/autotest/t003_test.py
index 70b4ed9928..e547ad2147 100644
--- a/autotest/t003_test.py
+++ b/autotest/t003_test.py
@@ -112,7 +112,7 @@ def test_load_nam_mt_nonexistant_file():
     test_loadoc_nstpfail()
     test_load_nam_mf_nonexistant_file()
     test_load_nam_mt_nonexistant_file()
-    # test_loadoc_lenfail()
-    # test_loadfreyberg()
-    # test_loadoahu()
-    # test_loadtwrip()
+    test_loadoc_lenfail()
+    test_loadfreyberg()
+    test_loadoahu()
+    test_loadtwrip()
diff --git a/autotest/t067_test_ulstrd.py b/autotest/t067_test_ulstrd.py
new file mode 100644
index 0000000000..7e771ea868
--- /dev/null
+++ b/autotest/t067_test_ulstrd.py
@@ -0,0 +1,137 @@
+import os
+import numpy as np
+import flopy
+
+tpth = os.path.join('temp', 't067')
+if not os.path.isdir(tpth):
+    os.makedirs(tpth)
+
+
+def test_ulstrd():
+
+    # Create an original model and then manually modify to use
+    # advanced list reader capabilities
+    ws = tpth
+    nlay = 1
+    nrow = 10
+    ncol = 10
+    nper = 3
+
+    # create the ghbs
+    ghbra = flopy.modflow.ModflowGhb.get_empty(20)
+    l = 0
+    for i in range(nrow):
+        ghbra[l] = (0, i, 0, 1., 100. + i)
+        l += 1
+        ghbra[l] = (0, i, ncol - 1, 1., 200. + i)
+        l += 1
+    ghbspd = {0: ghbra}
+
+    # create the drains
+    drnra = flopy.modflow.ModflowDrn.get_empty(2)
+    drnra[0] = (0, 1, int(ncol / 2), .5, 55.)
+    drnra[1] = (0, 2, int(ncol / 2), .5, 75.)
+    drnspd = {0: drnra}
+
+    # create the wells
+    welra = flopy.modflow.ModflowWel.get_empty(2)
+    welra[0] = (0, 1, 1, -5.)
+    welra[1] = (0, nrow - 3, ncol - 3, -10.)
+    welspd = {0: welra}
+
+    m = flopy.modflow.Modflow(modelname='original', model_ws=ws,
+                              exe_name='mf2005')
+    dis = flopy.modflow.ModflowDis(m, nlay=nlay, nrow=nrow, ncol=ncol,
+                                   nper=nper)
+    bas = flopy.modflow.ModflowBas(m)
+    lpf = flopy.modflow.ModflowLpf(m)
+    ghb = flopy.modflow.ModflowGhb(m, stress_period_data=ghbspd)
+    drn = flopy.modflow.ModflowDrn(m, stress_period_data=drnspd)
+    wel = flopy.modflow.ModflowWel(m, stress_period_data=welspd)
+    pcg = flopy.modflow.ModflowPcg(m)
+    oc = flopy.modflow.ModflowOc(m)
+    m.add_external('original.drn.dat', 71)
+    m.add_external('original.wel.bin', 72, binflag=True, output=False)
+    m.write_input()
+
+    # rewrite ghb
+    fname = os.path.join(ws, 'original.ghb')
+    with open(fname, 'w') as f:
+        f.write('{} {}\n'.format(ghbra.shape[0], 0))
+        for kper in range(nper):
+            f.write('{} {}\n'.format(ghbra.shape[0], 0))
+            f.write('open/close original.ghb.dat\n')
+
+    # write ghb list
+    sfacghb = 5
+    fname = os.path.join(ws, 'original.ghb.dat')
+    with open(fname, 'w') as f:
+        f.write('sfac {}\n'.format(sfacghb))
+        for k, i, j, stage, cond in ghbra:
+            f.write('{} {} {} {} {}\n'.format(k + 1, i + 1, j + 1, stage, cond))
+
+    # rewrite drn
+    fname = os.path.join(ws, 'original.drn')
+    with open(fname, 'w') as f:
+        f.write('{} {}\n'.format(drnra.shape[0], 0))
+        for kper in range(nper):
+            f.write('{} {}\n'.format(drnra.shape[0], 0))
+            f.write('external 71\n')
+
+    # write drn list
+    sfacdrn = 1.5
+    fname = os.path.join(ws, 'original.drn.dat')
+    with open(fname, 'w') as f:
+        for kper in range(nper):
+            f.write('sfac {}\n'.format(sfacdrn))
+            for k, i, j, stage, cond in drnra:
+                f.write(
+                    '{} {} {} {} {}\n'.format(k + 1, i + 1, j + 1, stage, cond))
+
+    # rewrite wel
+    fname = os.path.join(ws, 'original.wel')
+    with open(fname, 'w') as f:
+        f.write('{} {}\n'.format(drnra.shape[0], 0))
+        for kper in range(nper):
+            f.write('{} {}\n'.format(drnra.shape[0], 0))
+            f.write('external 72 (binary)\n')
+
+    # create the wells, but use an all float dtype to write a binary file
+    # use one-based values
+    weldt = np.dtype([('k', '<f4'), ('i', '<f4'), ('j', '<f4'), ('q', '<f4'), ])
+    welra = np.recarray(2, dtype=weldt)
+    welra[0] = (1, 2, 2, -5.)
+    welra[1] = (1, nrow - 2, ncol - 2, -10.)
+    fname = os.path.join(ws, 'original.wel.bin')
+    with open(fname, 'wb') as f:
+        welra.tofile(f)
+        welra.tofile(f)
+        welra.tofile(f)
+
+    # no need to run the model
+    #success, buff = m.run_model(silent=True)
+    #assert success, 'model did not terminate successfully'
+
+    # the m2 model will load all of these external files, possibly using sfac
+    # and just create regular list input files for wel, drn, and ghb
+    fname = 'original.nam'
+    m2 = flopy.modflow.Modflow.load(fname, model_ws=ws, verbose=False)
+    m2.name = 'new'
+    m2.write_input()
+
+    originalghbra = m.ghb.stress_period_data[0].copy()
+    originalghbra['cond'] *= sfacghb
+    assert np.array_equal(originalghbra, m2.ghb.stress_period_data[0])
+
+    originaldrnra = m.drn.stress_period_data[0].copy()
+    originaldrnra['cond'] *= sfacdrn
+    assert np.array_equal(originaldrnra, m2.drn.stress_period_data[0])
+
+    originalwelra = m.wel.stress_period_data[0].copy()
+    assert np.array_equal(originalwelra, m2.wel.stress_period_data[0])
+
+    return
+
+
+if __name__ == '__main__':
+    test_ulstrd()
diff --git a/examples/data/parameters/Oahu_02.lpf b/examples/data/parameters/Oahu_02.lpf
index 20732b7cac..a7c6dcdc7a 100644
--- a/examples/data/parameters/Oahu_02.lpf
+++ b/examples/data/parameters/Oahu_02.lpf
@@ -9,3 +9,6 @@ LAY1_HK  HK  1.0  1
 1 HKL1 ALL
 LAY1_VKA VK  1.0  1
 1 VKAL1 ALL
+1
+1
+1
diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py
index cc609eb5d0..7450b7bb3a 100644
--- a/flopy/modflow/mf.py
+++ b/flopy/modflow/mf.py
@@ -429,14 +429,17 @@ def write_name_file(self):
         f_nam.write('{}'.format(self.get_name_file_entries()))
 
         # write the external files
-        for u, f, b in zip(self.external_units, self.external_fnames,
-                           self.external_binflag):
+        for u, f, b, o in zip(self.external_units, self.external_fnames,
+                              self.external_binflag, self.external_output):
             if u == 0:
                 continue
-            # fr = os.path.relpath(f, self.model_ws)
+            replace_text = ''
+            if o:
+                replace_text = 'REPLACE'
             if b:
                 f_nam.write(
-                    'DATA(BINARY)   {0:5d}  '.format(u) + f + ' REPLACE\n')
+                    'DATA(BINARY)   {0:5d}  '.format(u) + f +
+                     replace_text + '\n')
             else:
                 f_nam.write('DATA           {0:5d}  '.format(u) + f + '\n')
 
diff --git a/flopy/modflow/mfchd.py b/flopy/modflow/mfchd.py
index eed8f9e411..894acec57f 100644
--- a/flopy/modflow/mfchd.py
+++ b/flopy/modflow/mfchd.py
@@ -192,6 +192,10 @@ def get_default_dtype(structured=True):
                               ("ehead", np.float32)])
         return dtype
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['shead', 'ehead']
+
     @staticmethod
     def load(f, model, nper=None, ext_unit_dict=None):
         """
diff --git a/flopy/modflow/mfdrn.py b/flopy/modflow/mfdrn.py
index d800c85890..f75715b999 100644
--- a/flopy/modflow/mfdrn.py
+++ b/flopy/modflow/mfdrn.py
@@ -241,6 +241,10 @@ def get_empty(ncells=0, aux_names=None, structured=True, is_drt=False):
             dtype = Package.add_to_dtype(dtype, aux_names, np.float32)
         return create_empty_recarray(ncells, dtype, default_value=-1.0E+10)
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['cond']
+
     @staticmethod
     def load(f, model, nper=None, ext_unit_dict=None, check=True):
         """
diff --git a/flopy/modflow/mfghb.py b/flopy/modflow/mfghb.py
index a3e309abda..6bb55bbe19 100644
--- a/flopy/modflow/mfghb.py
+++ b/flopy/modflow/mfghb.py
@@ -226,6 +226,10 @@ def get_default_dtype(structured=True):
                               ("cond", np.float32)])
         return dtype
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['cond']
+
     @staticmethod
     def load(f, model, nper=None, ext_unit_dict=None, check=True):
         """
diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py
index cfc7b7d5e5..3c694fb865 100644
--- a/flopy/modflow/mfhfb.py
+++ b/flopy/modflow/mfhfb.py
@@ -225,6 +225,10 @@ def get_default_dtype(structured=True):
             assert not structured, 'is there an unstructured HFB???'
         return dtype
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['hydchr']
+
     @staticmethod
     def load(f, model, ext_unit_dict=None):
         """
@@ -294,7 +298,9 @@ def load(f, model, ext_unit_dict=None):
         # data set 2 and 3
         if nphfb > 0:
             dt = ModflowHfb.get_empty(1).dtype
-            pak_parms = mfparbc.load(f, nphfb, dt, model.verbose)
+            pak_parms = mfparbc.load(f, nphfb, dt, model,
+                                     ext_unit_dict=ext_unit_dict,
+                                     verbose=model.verbose)
         # data set 4
         bnd_output = None
         if nhfbnp > 0:
@@ -328,9 +334,6 @@ def load(f, model, ext_unit_dict=None):
                 iname = 'static'
                 par_dict, current_dict = pak_parms.get(pname)
                 data_dict = current_dict[iname]
-                # print par_dict
-                # print data_dict
-
                 par_current = ModflowHfb.get_empty(par_dict['nlst'])
 
                 #
@@ -344,6 +347,7 @@ def load(f, model, ext_unit_dict=None):
 
                 # fill current parameter data (par_current)
                 for ibnd, t in enumerate(data_dict):
+                    t = tuple(t)
                     par_current[ibnd] = tuple(t[:len(par_current.dtype.names)])
 
                 # convert indices to zero-based
diff --git a/flopy/modflow/mfparbc.py b/flopy/modflow/mfparbc.py
index f208a8f477..e22b90bddf 100644
--- a/flopy/modflow/mfparbc.py
+++ b/flopy/modflow/mfparbc.py
@@ -5,7 +5,7 @@
 """
 
 import numpy as np
-from ..utils.flopy_io import line_strip
+from ..utils.flopy_io import line_strip, ulstrd
 
 
 class ModflowParBc(object):
@@ -42,7 +42,7 @@ def get(self, fkey):
         return None
 
     @staticmethod
-    def load(f, npar, dt, verbose=False):
+    def load(f, npar, dt, model, ext_unit_dict=None, verbose=False):
         """
         Load bc property parameters from an existing bc package
         that uses list data (e.g. WEL, RIV, etc.).
@@ -101,22 +101,12 @@ def load(f, npar, dt, verbose=False):
                         instnam = t[0].lower()
                     else:
                         instnam = 'static'
-                    bcinst = []
-                    for nw in range(nlst):
-                        line = f.readline()
-                        t = line_strip(line).split()
-                        bnd = []
-                        for jdx in range(nitems):
-                            # if jdx < 3:
-                            if issubclass(dt[jdx].type, np.integer):
-                                # conversion to zero-based occurs in package load method in mbase.
-                                bnd.append(int(t[jdx]))
-                            else:
-                                try:
-                                    bnd.append(float(t[jdx]))
-                                except IndexError:
-                                    bnd.append(1.)
-                        bcinst.append(bnd)
+
+                    ra = np.zeros(nlst, dtype=dt)
+                    #todo: if sfac is used for parameter definition, then
+                    # the empty list on the next line needs to be the package
+                    # get_sfac_columns
+                    bcinst = ulstrd(f, nlst, ra, model, [], ext_unit_dict)
                     pinst[instnam] = bcinst
                 bc_parms[parnam] = [{'partyp': partyp, 'parval': parval,
                                      'nlst': nlst, 'timevarying': timeVarying},
diff --git a/flopy/modflow/mfriv.py b/flopy/modflow/mfriv.py
index ea54173728..090f093ba5 100644
--- a/flopy/modflow/mfriv.py
+++ b/flopy/modflow/mfriv.py
@@ -253,6 +253,10 @@ def get_default_dtype(structured=True):
 
         return dtype
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['cond']
+
     def ncells(self):
         # Return the  maximum number of cells that have river
         # (developed for MT3DMS SSM package)
diff --git a/flopy/modflow/mfstr.py b/flopy/modflow/mfstr.py
index 555352973b..c47010ce37 100644
--- a/flopy/modflow/mfstr.py
+++ b/flopy/modflow/mfstr.py
@@ -619,7 +619,8 @@ def load(f, model, nper=None, ext_unit_dict=None):
         # read parameter data
         if npstr > 0:
             dt = ModflowStr.get_empty(1, aux_names=aux_names).dtype
-            pak_parms = mfparbc.load(f, npstr, dt, model.verbose)
+            pak_parms = mfparbc.load(f, npstr, dt, model, ext_unit_dict,
+                                     model.verbose)
 
         if nper is None:
             nrow, ncol, nlay, nper = model.get_nrow_ncol_nlay_nper()
diff --git a/flopy/modflow/mfwel.py b/flopy/modflow/mfwel.py
index 40b528c02c..a1b44fc23a 100644
--- a/flopy/modflow/mfwel.py
+++ b/flopy/modflow/mfwel.py
@@ -331,6 +331,10 @@ def get_empty(ncells=0, aux_names=None, structured=True):
             dtype = Package.add_to_dtype(dtype, aux_names, np.float32)
         return create_empty_recarray(ncells, dtype, default_value=-1.0E+10)
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['flux']
+
     @staticmethod
     def load(f, model, nper=None, ext_unit_dict=None, check=True):
         """
diff --git a/flopy/pakbase.py b/flopy/pakbase.py
index b15b3b79a7..e6f1d36e29 100644
--- a/flopy/pakbase.py
+++ b/flopy/pakbase.py
@@ -19,6 +19,7 @@
 from .modflow.mfparbc import ModflowParBc as mfparbc
 from .utils import Util2d, Util3d, Transient2d, MfList, check
 from .utils import OptionBlock
+from .utils.flopy_io import ulstrd
 
 
 class PackageInterface(object):
@@ -307,6 +308,15 @@ def add_to_dtype(dtype, field_names, field_types):
             newdtypes.append((str(field_name), field_type))
         return np.dtype(newdtypes)
 
+    @staticmethod
+    def get_sfac_columns():
+        """
+        This should be overriden for individual packages that support an
+        sfac multiplier for individual list columns
+
+        """
+        return []
+
     def check(self, f=None, verbose=True, level=1):
         """
         Check package data for common errors.
@@ -790,12 +800,15 @@ def load(f, model, pak_type, ext_unit_dict=None, **kwargs):
         elif 'flopy.modflow.mfchd.modflowchd'.lower() in pak_type_str:
             partype = ['shead', 'ehead']
 
+        # get the list columns that should be scaled with sfac
+        sfac_columns = pak_type.get_sfac_columns()
+
         # read parameter data
         if nppak > 0:
             dt = pak_type.get_empty(1, aux_names=aux_names,
                                     structured=model.structured).dtype
-            pak_parms = mfparbc.load(f, nppak, dt, model.verbose)
-            # pak_parms = mfparbc.load(f, nppak, len(dt.names))
+            pak_parms = mfparbc.load(f, nppak, dt, model, ext_unit_dict,
+                                     model.verbose)
 
         if nper is None:
             nrow, ncol, nlay, nper = model.get_nrow_ncol_nlay_nper()
@@ -827,69 +840,8 @@ def load(f, model, pak_type, ext_unit_dict=None, **kwargs):
             elif itmp > 0:
                 current = pak_type.get_empty(itmp, aux_names=aux_names,
                                              structured=model.structured)
-                for ibnd in range(itmp):
-                    line = f.readline()
-                    if "open/close" in line.lower():
-                        binary = False
-                        if '(binary)' in line.lower():
-                            binary = True
-                        # need to strip out existing path seps and
-                        # replace current-system path seps
-                        raw = line.strip().split()
-                        fname = raw[1]
-                        if '/' in fname:
-                            raw = fname.split('/')
-                        elif '\\' in fname:
-                            raw = fname.split('\\')
-                        else:
-                            raw = [fname]
-                        fname = os.path.join(*raw)
-                        oc_filename = os.path.join(model.model_ws, fname)
-                        msg = 'Package.load() error: open/close filename ' + \
-                              oc_filename + ' not found'
-                        assert os.path.exists(oc_filename), msg
-                        try:
-                            if binary:
-                                dtype2 = []
-                                for name in current.dtype.names:
-                                    dtype2.append((name, np.float32))
-                                dtype2 = np.dtype(dtype2)
-                                d = np.fromfile(oc_filename,
-                                                dtype=dtype2,
-                                                count=itmp)
-                                current = np.array(d, dtype=current.dtype)
-                            else:
-                                cd = current.dtype
-                                current = np.loadtxt(oc_filename).transpose()
-                                if current.ndim == 1:
-                                    current = np.atleast_2d(
-                                        current).transpose()
-                                current = np.core.records.fromarrays(current,
-                                                                     dtype=cd)
-                            current = current.view(np.recarray)
-                        except Exception as e:
-                            msg = 'Package.load() error loading ' + \
-                                  'open/close file ' + oc_filename + \
-                                  ': ' + str(e)
-                            raise Exception(msg)
-                        msg = 'Package.load() error: open/close ' + \
-                              'recarray from file ' + oc_filename + \
-                              ' shape (' + str(current.shape) + \
-                              ') does not match itmp: {:d}'.format(itmp)
-                        assert current.shape[0] == itmp, msg
-                        break
-                    try:
-                        t = line.strip().split()
-                        current[ibnd] = tuple(t[:len(current.dtype.names)])
-                    except:
-                        t = []
-                        for ivar in range(len(current.dtype.names)):
-                            istart = ivar * 10
-                            istop = istart + 10
-                            t.append(line[istart:istop])
-                        current[ibnd] = tuple(t[:len(current.dtype.names)])
-
-                # convert indices to zero-based
+                current = ulstrd(f, itmp, current, model, sfac_columns,
+                                 ext_unit_dict)
                 if model.structured:
                     current['k'] -= 1
                     current['i'] -= 1
@@ -935,6 +887,7 @@ def load(f, model, pak_type, ext_unit_dict=None, **kwargs):
 
                 # fill current parameter data (par_current)
                 for ibnd, t in enumerate(data_dict):
+                    t = tuple(t)
                     par_current[ibnd] = tuple(t[:len(par_current.dtype.names)])
 
                 if model.structured:
diff --git a/flopy/utils/flopy_io.py b/flopy/utils/flopy_io.py
index 3c33c1844d..531abc7ddf 100755
--- a/flopy/utils/flopy_io.py
+++ b/flopy/utils/flopy_io.py
@@ -1,6 +1,7 @@
 """
 Module for input/output utilities
 """
+import os
 import sys
 import numpy as np
 
@@ -323,3 +324,134 @@ def get_url_text(url, error_msg=None):
         if error_msg is not None:
             print(error_msg)
         return
+
+
+def ulstrd(f, nlist, ra, model, sfac_columns, ext_unit_dict):
+    """
+    Read a list and allow for open/close, binary, external, sfac, etc.
+
+    Parameters
+    ----------
+    f : file handle
+        file handle for where the list is being read from
+    nlist : int
+        size of the list (number of rows) to read
+    ra : np.recarray
+        A record array of the correct size that will be filled with the list
+    model : model object
+        The model object (of type :class:`flopy.modflow.mf.Modflow`) to
+        which this package will be added.
+    sfac_columns : list
+        A list of strings containing the column names to scale by sfac
+    ext_unit_dict : dictionary, optional
+        If the list in the file is specified using EXTERNAL,
+        then in this case ext_unit_dict is required, which can be
+        constructed using the function
+        :class:`flopy.utils.mfreadnam.parsenamefile`.
+
+    Returns
+    -------
+
+    """
+
+    # initialize variables
+    line = f.readline()
+    sfac = 1.
+    binary = False
+    ncol = len(ra.dtype.names)
+    line_list = line.strip().split()
+    close_the_file = False
+    file_handle = f
+    mode = 'r'
+
+    # check for external
+    if line.strip().lower().startswith('external'):
+        inunit = int(line_list[1])
+        errmsg = 'Could not find a file for unit {}'.format(inunit)
+        if ext_unit_dict is not None:
+            if inunit in ext_unit_dict:
+                namdata = ext_unit_dict[inunit]
+                file_handle = namdata.filehandle
+            else:
+                raise IOError(errmsg)
+        else:
+            raise IOError(errmsg)
+        if namdata.filetype == 'DATA(BINARY)':
+            binary = True
+        if not binary:
+            line = file_handle.readline()
+
+    # or check for open/close
+    elif line.strip().lower().startswith('open/close'):
+        raw = line.strip().split()
+        fname = raw[1]
+        if '/' in fname:
+            raw = fname.split('/')
+        elif '\\' in fname:
+            raw = fname.split('\\')
+        else:
+            raw = [fname]
+        fname = os.path.join(*raw)
+        oc_filename = os.path.join(model.model_ws, fname)
+        msg = 'Package.load() error: open/close filename ' + \
+              oc_filename + ' not found'
+        assert os.path.exists(oc_filename), msg
+        if '(binary)' in line.lower():
+            binary = True
+            mode = 'rb'
+        file_handle = open(oc_filename, mode)
+        close_the_file = True
+        if not binary:
+            line = file_handle.readline()
+
+    # check for scaling factor
+    if not binary:
+        line_list = line.strip().split()
+        if line.strip().lower().startswith('sfac'):
+            sfac = float(line_list[1])
+            line = file_handle.readline()
+
+    # fast binary read fromfile
+    if binary:
+        dtype2 = []
+        for name in ra.dtype.names:
+            dtype2.append((name, np.float32))
+        dtype2 = np.dtype(dtype2)
+        d = np.fromfile(file_handle, dtype=dtype2, count=nlist)
+        ra = np.array(d, dtype=ra.dtype)
+        ra = ra.view(np.recarray)
+
+    # else, read ascii
+    else:
+
+        for ii in range(nlist):
+
+            # first line was already read
+            if ii != 0:
+                line = file_handle.readline()
+
+            if model.free_format_input:
+                # whitespace separated
+                t = line.strip().split()
+                if len(t) < ncol:
+                    t = t + (ncol - len(t)) * [0.0]
+                else:
+                    t = t[:ncol]
+                t = tuple(t)
+                ra[ii] = t
+            else:
+                # fixed format
+                t = read_fixed_var(line, ncol=ncol)
+                t = tuple(t)
+                ra[ii] = t
+
+    # scale the data and check
+    for column_name in sfac_columns:
+        ra[column_name] *= sfac
+        if 'auxsfac' in ra.dtype.names:
+            ra[column_name] *= ra['auxsfac']
+
+    if close_the_file:
+        file_handle.close()
+
+    return ra