diff --git a/autotest/t007_test.py b/autotest/t007_test.py
index 40e5537551..a81ce61682 100644
--- a/autotest/t007_test.py
+++ b/autotest/t007_test.py
@@ -163,6 +163,19 @@ def export_shapefile_modelgrid_override(namfile):
         raise Exception(msg)
 
 
+def test_output_helper_shapefile_export():
+    ws = os.path.join('..', 'examples', 'data', 'freyberg_multilayer_transient')
+    name = 'freyberg.nam'
+
+    ml = flopy.modflow.Modflow.load(name, model_ws=ws)
+
+    head = flopy.utils.HeadFile(os.path.join(ws, 'freyberg.hds'))
+    cbc = flopy.utils.CellBudgetFile(os.path.join(ws, "freyberg.cbc"))
+    flopy.export.utils.output_helper(os.path.join('temp', 'test.shp'), ml,
+                                     {'HDS': head, 'cbc': cbc},
+                                     mflay=1, kper=10)
+
+
 def test_freyberg_export():
     from flopy.discretization import StructuredGrid
     namfile = 'freyberg.nam'
@@ -1119,6 +1132,174 @@ def check_vertices():
     plt.close()
 
 
+def test_mapview_plot_bc():
+    from matplotlib.collections import QuadMesh, PatchCollection
+    import matplotlib.pyplot as plt
+
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join("..", "examples", "data", "mf6",
+                            "test003_gwfs_disv")
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+    ml6 = sim.get_model("gwf_1")
+    ml6.modelgrid.set_coord_info(angrot=-14)
+    mapview = flopy.plot.PlotMapView(model=ml6)
+    mapview.plot_bc('CHD')
+    ax = mapview.ax
+
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, PatchCollection):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6', 'test045_lake2tr')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+
+    ml6 = sim.get_model("lakeex2a")
+    mapview = flopy.plot.PlotMapView(model=ml6)
+    mapview.plot_bc('LAK')
+    mapview.plot_bc("SFR")
+
+    ax = mapview.ax
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, QuadMesh):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6',
+                            'test006_2models_mvr')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+
+    ml6 = sim.get_model("parent")
+    ml6c = sim.get_model('child')
+    ml6c.modelgrid.set_coord_info(xoff=700, yoff=0, angrot=0)
+
+    mapview = flopy.plot.PlotMapView(model=ml6)
+    mapview.plot_bc("MAW")
+
+    mapview2 = flopy.plot.PlotMapView(model=ml6c, ax=mapview.ax)
+    mapview2.plot_bc("MAW")
+    ax = mapview2.ax
+
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, QuadMesh):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6',
+                            'test001e_UZF_3lay')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+    ml6 = sim.get_model("gwf_1")
+
+    mapview = flopy.plot.PlotMapView(model=ml6)
+    mapview.plot_bc("UZF")
+
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, QuadMesh):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+
+def test_crosssection_plot_bc():
+    from matplotlib.collections import PatchCollection
+    import matplotlib.pyplot as plt
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join("..", "examples", "data", "mf6",
+                            "test003_gwfs_disv")
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+    ml6 = sim.get_model("gwf_1")
+    xc = flopy.plot.PlotCrossSection(ml6, line={'line': ([0, 5.5],
+                                                         [10, 5.5])})
+    xc.plot_bc('CHD')
+    ax = xc.ax
+
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, PatchCollection):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6', 'test045_lake2tr')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+
+    ml6 = sim.get_model("lakeex2a")
+    xc = flopy.plot.PlotCrossSection(ml6, line={'row': 10})
+    xc.plot_bc('LAK')
+    xc.plot_bc("SFR")
+
+    ax = xc.ax
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, PatchCollection):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6',
+                            'test006_2models_mvr')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+
+    ml6 = sim.get_model("parent")
+    xc = flopy.plot.PlotCrossSection(ml6, line={'column': 1})
+    xc.plot_bc("MAW")
+
+    ax = xc.ax
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, PatchCollection):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6',
+                            'test001e_UZF_3lay')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+    ml6 = sim.get_model("gwf_1")
+
+    xc = flopy.plot.PlotCrossSection(ml6, line={"row": 0})
+    xc.plot_bc("UZF")
+
+    ax = xc.ax
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, PatchCollection):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+
 def test_tricontour_NaN():
     from flopy.plot import PlotMapView
     import numpy as np
@@ -1258,43 +1439,6 @@ def test_netcdf_classmethods():
     assert len(diff) == 0, str(diff)
 
 
-# def test_netcdf_overloads():
-#     import os
-#     import flopy
-#     nam_file = "freyberg.nam"
-#     model_ws = os.path.join('..', 'examples', 'data', 'freyberg_multilayer_transient')
-#     ml = flopy.modflow.Modflow.load(nam_file,model_ws=model_ws,check=False,
-#                                     verbose=False,load_only=[])
-#
-#     f = ml.export(os.path.join("temp","freyberg.nc"))
-#     fzero = flopy.export.NetCdf.zeros_like(f)
-#     assert fzero.nc.variables["model_top"][:].sum() == 0
-#     print(f.nc.variables["model_top"][0,:])
-#     fplus1 = f + 1
-#     assert fplus1.nc.variables["model_top"][0,0] == f.nc.variables["model_top"][0,0] + 1
-#     assert (f + fplus1).nc.variables["model_top"][0,0] ==\
-#            f.nc.variables["model_top"][0,0] + \
-#            fplus1.nc.variables["model_top"][0,0]
-#
-#     fminus1 = f - 1
-#     assert fminus1.nc.variables["model_top"][0,0] == f.nc.variables["model_top"][0,0] - 1
-#     assert (f - fminus1).nc.variables["model_top"][0,0]==\
-#            f.nc.variables["model_top"][0,0] - \
-#            fminus1.nc.variables["model_top"][0,0]
-#
-#     ftimes2 = f * 2
-#     assert ftimes2.nc.variables["model_top"][0,0] == f.nc.variables["model_top"][0,0] * 2
-#     assert (f * ftimes2).nc.variables["model_top"][0,0] ==\
-#             f.nc.variables["model_top"][0,0] * \
-#            ftimes2.nc.variables["model_top"][0,0]
-#
-#     fdiv2 = f / 2
-#     assert fdiv2.nc.variables["model_top"][0,0] == f.nc.variables["model_top"][0,0] / 2
-#     assert (f / fdiv2).nc.variables["model_top"][0,0] == \
-#          f.nc.variables["model_top"][0,0] / \
-#            fdiv2.nc.variables["model_top"][0,0]
-#
-#     assert f.nc.variables["ibound"][0,0,0] == 1
 def test_wkt_parse():
     """Test parsing of Coordinate Reference System parameters
     from well-known-text in .prj files."""
@@ -1458,15 +1602,14 @@ def test_export_contourf():
 def main():
     # test_shapefile()
     # test_shapefile_ibound()
+    # test_netcdf_classmethods()
 
-    test_netcdf_classmethods()
-
-    for namfile in namfiles:
-        export_mf2005_netcdf(namfile)
-        export_shapefile(namfile)
+    # for namfile in namfiles:
+    #    export_mf2005_netcdf(namfile)
+    #    export_shapefile(namfile)
 
-    for namfile in namfiles[0:2]:
-        export_shapefile_modelgrid_override(namfile)
+    # for namfile in namfiles[0:2]:
+    #     export_shapefile_modelgrid_override(namfile)
 
     # test_netcdf_overloads()
     # test_netcdf_classmethods()
@@ -1480,7 +1623,7 @@ def main():
     # test_vertex_model_dot_plot()
     # test_sr_with_Map()
     # test_modelgrid_with_PlotMapView()
-    test_epsgs()
+    # test_epsgs()
     # test_sr_scaling()
     # test_read_usgs_model_reference()
     # test_dynamic_xll_yll()
@@ -1500,7 +1643,9 @@ def main():
     # test_export_contourf()
     # test_sr()
     # test_shapefile_polygon_closed()
-
+    test_mapview_plot_bc()
+    test_crosssection_plot_bc()
+    test_output_helper_shapefile_export()
 
 if __name__ == '__main__':
 
diff --git a/autotest/t012_test.py b/autotest/t012_test.py
index e0d65a9a6d..70e71026dd 100644
--- a/autotest/t012_test.py
+++ b/autotest/t012_test.py
@@ -61,7 +61,7 @@ def test_mf2000_p07():
     pth = os.path.join(pth2000, 'P07')
     namfile = 'p7mf2k.nam'
     mf = flopy.modflow.Modflow.load(namfile, model_ws=pth,
-                                    version='mf2k', verbose=True,
+                                    verbose=True,
                                     exe_name=mf2k_exe)
 
     cpth = os.path.join(newpth, 'P07_2K')
diff --git a/flopy/export/utils.py b/flopy/export/utils.py
index b212b0630c..4b4a9eba89 100644
--- a/flopy/export/utils.py
+++ b/flopy/export/utils.py
@@ -278,6 +278,10 @@ def output_helper(f, ml, oudic, **kwargs):
         modelgrid : flopy.discretizaiton.Grid
             user supplied model grid instance that will be used for export
             in lieu of the models model grid instance
+        mflay : int
+            zero based model layer which can be used in shapefile exporting
+        kper : int
+            zero based stress period which can be used for shapefile exporting
 
     Returns
     -------
@@ -294,6 +298,8 @@ def output_helper(f, ml, oudic, **kwargs):
     forgive = kwargs.pop("forgive", False)
     kwargs.pop("suffix", None)
     mask_vals = []
+    mflay = kwargs.pop('mflay', None)
+    kper = kwargs.pop('kper', None)
     if "masked_vals" in kwargs:
         mask_vals = kwargs.pop("masked_vals")
     if len(kwargs) > 0 and logger is not None:
@@ -427,11 +433,65 @@ def output_helper(f, ml, oudic, **kwargs):
                                                    zonebud.flux,
                                                    logger)
 
-            # todo: write the zone array to standard output
+            # write the zone array to standard output
             _add_output_nc_variable(f, times, shape3d, zonebud,
                                     "budget_zones", logger=logger,
                                     mask_vals=mask_vals,
                                     mask_array3d=mask_array3d)
+
+    elif isinstance(f, str) and f.endswith('.shp'):
+        attrib_dict = {}
+        for _, out_obj in oudic.items():
+
+            if isinstance(out_obj, HeadFile) or \
+                    isinstance(out_obj, FormattedHeadFile) or \
+                    isinstance(out_obj, UcnFile):
+                if isinstance(out_obj, UcnFile):
+                    attrib_name = 'conc'
+                else:
+                    attrib_name = 'head'
+                plotarray = np.atleast_3d(out_obj.get_alldata()
+                                          .transpose()).transpose()
+
+
+                for per in range(plotarray.shape[0]):
+                    for k in range(plotarray.shape[1]):
+                        if kper is not None and per != kper:
+                            continue
+                        if mflay is not None and k != mflay:
+                            continue
+                        name = attrib_name + '{}_{}'.format(per, k)
+                        attrib_dict[name] = plotarray[per][k]
+
+            elif isinstance(out_obj, CellBudgetFile):
+                names = out_obj.get_unique_record_names(decode=True)
+
+                for attrib_name in names:
+                    plotarray = np.atleast_3d(out_obj.get_data(
+                                              text=attrib_name,
+                                              full3D=True))
+
+                    attrib_name = attrib_name.strip()
+                    if attrib_name == "FLOW RIGHT FACE":
+                        attrib_name = 'FRF'
+                    elif attrib_name == "FLOW FRONT FACE":
+                        attrib_name = "FFF"
+                    elif attrib_name == "FLOW LOWER FACE":
+                        attrib_name = "FLF"
+                    else:
+                        pass
+                    for per in range(plotarray.shape[0]):
+                        for k in range(plotarray.shape[1]):
+                            if kper is not None and per != kper:
+                                continue
+                            if mflay is not None and k != mflay:
+                                continue
+                            name = attrib_name + '{}_{}'.format(per, k)
+                            attrib_dict[name] = plotarray[per][k]
+
+        if attrib_dict:
+            shapefile_utils.write_grid_shapefile(f, ml.modelgrid, attrib_dict)
+
     else:
         if logger:
             logger.lraise("unrecognized export argument:{0}".format(f))
diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py
index 044017189b..4e84ef8ddf 100644
--- a/flopy/modflow/mf.py
+++ b/flopy/modflow/mf.py
@@ -682,6 +682,8 @@ def load(f, version='mf2005', exe_name='mf2005.exe', verbose=False,
         if 'NWT' in ext_pkg_d or 'UPW' in ext_pkg_d:
             version = 'mfnwt'
         if 'GLOBAL' in ext_pkg_d:
+            if version != "mf2k":
+                ml.glo = ModflowGlobal(ml)
             version = 'mf2k'
         if 'SMS' in ext_pkg_d:
             version = 'mfusg'
diff --git a/flopy/pakbase.py b/flopy/pakbase.py
index befc2b36a7..72b6bdadff 100644
--- a/flopy/pakbase.py
+++ b/flopy/pakbase.py
@@ -502,8 +502,8 @@ def parent(self, parent):
 
     @property
     def package_type(self):
-        if len(self.extension) > 0:
-            return self.extension[0]
+        if len(self.name) > 0:
+            return self.name[0].lower()
 
     @property
     def plotable(self):
diff --git a/flopy/plot/map.py b/flopy/plot/map.py
index 1f39d11c20..89cb3f7917 100644
--- a/flopy/plot/map.py
+++ b/flopy/plot/map.py
@@ -452,7 +452,7 @@ def plot_grid(self, **kwargs):
 
         return lc
 
-    def plot_bc(self, ftype=None, package=None, kper=0, color=None,
+    def plot_bc(self, name=None, package=None, kper=0, color=None,
                 plotAll=False, **kwargs):
         """
         Plot boundary conditions locations for a specific boundary
@@ -460,7 +460,7 @@ def plot_bc(self, ftype=None, package=None, kper=0, color=None,
 
         Parameters
         ----------
-        ftype : string
+        name : string
             Package name string ('WEL', 'GHB', etc.). (Default is None)
         package : flopy.modflow.Modflow package class instance
             flopy package class instance. (Default is None)
@@ -480,41 +480,66 @@ def plot_bc(self, ftype=None, package=None, kper=0, color=None,
         quadmesh : matplotlib.collections.QuadMesh
 
         """
+        if 'ftype' in kwargs and name is None:
+            name = kwargs.pop('ftype')
+
         # Find package to plot
         if package is not None:
             p = package
-            ftype = p.name[0]
+            name = p.name[0]
 
         elif self.model is not None:
-            if ftype is None:
+            if name is None:
                 raise Exception('ftype not specified')
-            ftype = ftype.upper()
-            p = self.model.get_package(ftype)
+            name = name.upper()
+            p = self.model.get_package(name)
 
         else:
             raise Exception('Cannot find package to plot')
 
         # trap for mf6 'cellid' vs mf2005 'k', 'i', 'j' convention
-        if p.parent.version == "mf6":
-            try:
-                mflist = p.stress_period_data.array[kper]
-            except Exception as e:
-                raise Exception("Not a list-style boundary package: " + str(e))
-            if mflist is None:
-                return
-            idx = np.array([list(i) for i in mflist['cellid']], dtype=int).T
+        if isinstance(p, list) or p.parent.version == "mf6":
+            if not isinstance(p, list):
+                p = [p]
+
+            idx = np.array([])
+            for pp in p:
+                if pp.package_type in ('lak', 'sfr', 'maw', 'uzf'):
+                    t = plotutil.advanced_package_bc_helper(pp, self.mg,
+                                                            kper)
+                else:
+                    try:
+                        mflist = pp.stress_period_data.array[kper]
+                    except Exception as e:
+                        raise Exception("Not a list-style boundary package: "
+                                        + str(e))
+                    if mflist is None:
+                        return
+
+                    t = np.array([list(i) for i in mflist['cellid']],
+                                     dtype=int).T
+
+                if len(idx) == 0:
+                    idx = np.copy(t)
+                else:
+                    idx = np.append(idx, t, axis=1)
+
         else:
             # modflow-2005 structured and unstructured grid
-            try:
-                mflist = p.stress_period_data[kper]
-            except Exception as e:
-                raise Exception("Not a list-style boundary package: " + str(e))
-            if mflist is None:
-                return
-            if len(self.mg.shape) == 3:
-                idx = [mflist['k'], mflist['i'], mflist['j']]
+            if p.package_type in ('uzf', 'lak'):
+                idx = plotutil.advanced_package_bc_helper(p, self.mg, kper)
             else:
-                idx = mflist['node']
+                try:
+                    mflist = p.stress_period_data[kper]
+                except Exception as e:
+                    raise Exception("Not a list-style boundary package: "
+                                    + str(e))
+                if mflist is None:
+                    return
+                if len(self.mg.shape) == 3:
+                    idx = [mflist['k'], mflist['i'], mflist['j']]
+                else:
+                    idx = mflist['node']
 
         nlay = self.mg.nlay
 
@@ -534,7 +559,7 @@ def plot_bc(self, ftype=None, package=None, kper=0, color=None,
         # set the colormap
         if color is None:
             # modflow 6 ftype fix, since multiple packages append _0, _1, etc:
-            key = ftype[:3].upper()
+            key = name[:3].upper()
             if key in plotutil.bc_color_dict:
                 c = plotutil.bc_color_dict[key]
             else:
diff --git a/flopy/plot/plotbase.py b/flopy/plot/plotbase.py
index 9636108d3f..ec54f5137f 100644
--- a/flopy/plot/plotbase.py
+++ b/flopy/plot/plotbase.py
@@ -304,7 +304,7 @@ def plot_grid(self, **kwargs):
 
         return col
 
-    def plot_bc(self, ftype=None, package=None, kper=0, color=None,
+    def plot_bc(self, name=None, package=None, kper=0, color=None,
                 head=None, **kwargs):
         """
         Plot boundary conditions locations for a specific boundary
@@ -312,7 +312,7 @@ def plot_bc(self, ftype=None, package=None, kper=0, color=None,
 
         Parameters
         ----------
-        ftype : string
+        name : string
             Package name string ('WEL', 'GHB', etc.). (Default is None)
         package : flopy.modflow.Modflow package class instance
             flopy package class instance. (Default is None)
@@ -334,51 +334,76 @@ def plot_bc(self, ftype=None, package=None, kper=0, color=None,
         patches : matplotlib.collections.PatchCollection
 
         """
+        if 'ftype' in kwargs and name is None:
+            name = kwargs.pop('ftype')
+
         # Find package to plot
         if package is not None:
             p = package
             ftype = p.name[0]
         elif self.model is not None:
-            if ftype is None:
+            if name is None:
                 raise Exception('ftype not specified')
-            ftype = ftype.upper()
-            p = self.model.get_package(ftype)
+            name = name.upper()
+            p = self.model.get_package(name)
         else:
             raise Exception('Cannot find package to plot')
 
         # trap for mf6 'cellid' vs mf2005 'k', 'i', 'j' convention
-        if p.parent.version == "mf6":
-            try:
-                mflist = p.stress_period_data.array[kper]
-            except Exception as e:
-                raise Exception("Not a list-style boundary package: " + str(e))
-            if mflist is None:
-                return
-            idx = np.array([list(i) for i in mflist['cellid']], dtype=int).T
+        if isinstance(p, list) or p.parent.version == "mf6":
+            if not isinstance(p, list):
+                p = [p]
+
+            idx = np.array([])
+            for pp in p:
+                if pp.package_type in ('lak', 'sfr', 'maw', 'uzf'):
+                    t = plotutil.advanced_package_bc_helper(pp, self.mg,
+                                                            kper)
+                else:
+                    try:
+                        mflist = pp.stress_period_data.array[kper]
+                    except Exception as e:
+                        raise Exception("Not a list-style boundary package: "
+                                        + str(e))
+                    if mflist is None:
+                        return
+
+                    t = np.array([list(i) for i in mflist['cellid']],
+                                 dtype=int).T
+
+                if len(idx) == 0:
+                    idx = np.copy(t)
+                else:
+                    idx = np.append(idx, t, axis=1)
+
         else:
             # modflow-2005 structured and unstructured grid
-            try:
-                mflist = p.stress_period_data[kper]
-            except Exception as e:
-                raise Exception("Not a list-style boundary package: " + str(e))
-            if mflist is None:
-                return
-            if len(self.mg.shape) == 3:
-                idx = [mflist['k'], mflist['i'], mflist['j']]
+            if p.package_type in ('uzf', 'lak'):
+                idx = plotutil.advanced_package_bc_helper(p, self.mg, kper)
             else:
-                idx = mflist['node']
+                try:
+                    mflist = p.stress_period_data[kper]
+                except Exception as e:
+                    raise Exception("Not a list-style boundary package: "
+                                    + str(e))
+                if mflist is None:
+                    return
+                if len(self.mg.shape) == 3:
+                    idx = [mflist['k'], mflist['i'], mflist['j']]
+                else:
+                    idx = mflist['node']
 
         # Plot the list locations, change this to self.mg.shape
-        if self.mg.grid_type == "vertex":
+        if len(self.mg.shape) != 3:
             plotarray = np.zeros((self.mg.nlay, self.mg.ncpl), dtype=np.int)
-            plotarray[idx] = 1
+            plotarray[tuple(idx)] = 1
         else:
             plotarray = np.zeros((self.mg.nlay, self.mg.nrow, self.mg.ncol), dtype=np.int)
             plotarray[idx[0], idx[1], idx[2]] = 1
 
         plotarray = np.ma.masked_equal(plotarray, 0)
         if color is None:
-            key = ftype[:3].upper()
+            key = name[:3].upper()
             if key in plotutil.bc_color_dict:
                 c = plotutil.bc_color_dict[key]
             else:
diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py
index aed55ee249..e75ca63a8b 100644
--- a/flopy/plot/plotutil.py
+++ b/flopy/plot/plotutil.py
@@ -288,8 +288,9 @@
 
 
 bc_color_dict = {'default': 'black', 'WEL': 'red', 'DRN': 'yellow',
-                 'RIV': 'green', 'GHB': 'cyan', 'CHD': 'navy',
-                 'STR': 'purple', 'SFR': 'blue'}
+                 'RIV': 'teal', 'GHB': 'cyan', 'CHD': 'navy',
+                 'STR': 'purple', 'SFR': 'teal', 'UZF': 'peru',
+                 'LAK': 'royalblue'}
 
 
 class PlotException(Exception):
@@ -2749,4 +2750,35 @@ def _depreciated_dis_handler(modelgrid, dis):
     return modelgrid
 
 
+def advanced_package_bc_helper(pkg, modelgrid, kper):
+    """
+    Helper function for plotting boundary conditions from "advanced" packages
 
+    Parameters
+    ----------
+    pkg : flopy Package objects
+    modelgrid : flopy.discretization.Grid object
+
+    Returns
+    -------
+    """
+    if pkg.package_type in ('sfr', 'uzf'):
+        if pkg.parent.version == 'mf6':
+            mflist = pkg.packagedata.array
+            idx = np.array([list(i) for i in mflist['cellid']], dtype=int).T
+        else:
+            iuzfbnd = pkg.iuzfbnd.array
+            idx = np.where(iuzfbnd != 0)
+            idx = np.append([[0] * idx[-1].size], idx, axis=0)
+    elif pkg.package_type in ('lak', 'maw'):
+        if pkg.parent.version == "mf6":
+            mflist = pkg.connectiondata.array
+            idx = np.array([list(i) for i in mflist['cellid']], dtype=int).T
+        else:
+            lakarr = pkg.lakarr.array[kper]
+            idx = np.where(lakarr != 0)
+            idx = np.array(idx)
+    else:
+        raise NotImplementedError("Pkg {} not implemented for bc plotting"
+                                  .format(pkg.package_type))
+    return idx