diff --git a/.docs/Notebooks/mfusg_transport_Ex1_1D.py b/.docs/Notebooks/mfusg_transport_tutorial01_1D.py
similarity index 98%
rename from .docs/Notebooks/mfusg_transport_Ex1_1D.py
rename to .docs/Notebooks/mfusg_transport_tutorial01_1D.py
index f289a66189..46cf29adcf 100644
--- a/.docs/Notebooks/mfusg_transport_Ex1_1D.py
+++ b/.docs/Notebooks/mfusg_transport_tutorial01_1D.py
@@ -41,8 +41,8 @@
 import numpy as np
 
 import flopy
-from flopy.mfusg import MfUsg, MfUsgBct, MfUsgLpf, MfUsgOc, MfUsgPcb, MfUsgSms
-from flopy.modflow import ModflowBas, ModflowDis
+from flopy.mfusg import MfUsg, MfUsgBas, MfUsgBct, MfUsgLpf, MfUsgOc, MfUsgPcb, MfUsgSms
+from flopy.modflow import ModflowDis
 from flopy.utils import HeadFile
 
 # -
@@ -95,7 +95,7 @@
 strt = np.full((nlay, nrow, ncol), 100.0)
 strt[:, :, 0] = 1100.0
 
-bas = ModflowBas(mf, ibound=ibound, strt=strt)
+bas = MfUsgBas(mf, ibound=ibound, strt=strt)
 # -
 
 # +
diff --git a/.docs/Notebooks/mfusg_transport_Ex2_Radial.py b/.docs/Notebooks/mfusg_transport_tutorial02_Radial.py
similarity index 96%
rename from .docs/Notebooks/mfusg_transport_Ex2_Radial.py
rename to .docs/Notebooks/mfusg_transport_tutorial02_Radial.py
index 22915173c7..92580d2cd8 100644
--- a/.docs/Notebooks/mfusg_transport_Ex2_Radial.py
+++ b/.docs/Notebooks/mfusg_transport_tutorial02_Radial.py
@@ -44,8 +44,8 @@
 import numpy as np
 
 import flopy
-from flopy.mfusg import MfUsg, MfUsgBct, MfUsgLpf, MfUsgOc, MfUsgSms, MfUsgWel
-from flopy.modflow import ModflowBas, ModflowChd, ModflowDis
+from flopy.mfusg import MfUsg, MfUsgBas, MfUsgBct, MfUsgLpf, MfUsgOc, MfUsgSms, MfUsgWel
+from flopy.modflow import ModflowChd, ModflowDis
 from flopy.utils import HeadFile
 
 # -
@@ -101,7 +101,7 @@
 strt[:, :, 0] = 20.0
 strt[:, :, -1] = 20.0
 
-bas = ModflowBas(mf, ibound=ibound, strt=strt)
+bas = MfUsgBas(mf, ibound=ibound, strt=strt)
 # -
 # +
 ipakcb = 50
diff --git a/.docs/Notebooks/mfusg_transport_Ex3_Conduit.py b/.docs/Notebooks/mfusg_transport_tutorial03_Conduit.py
similarity index 97%
rename from .docs/Notebooks/mfusg_transport_Ex3_Conduit.py
rename to .docs/Notebooks/mfusg_transport_tutorial03_Conduit.py
index c98221f375..3764b1999f 100644
--- a/.docs/Notebooks/mfusg_transport_Ex3_Conduit.py
+++ b/.docs/Notebooks/mfusg_transport_tutorial03_Conduit.py
@@ -38,8 +38,17 @@
 import numpy as np
 
 import flopy
-from flopy.mfusg import MfUsg, MfUsgBcf, MfUsgBct, MfUsgCln, MfUsgOc, MfUsgSms, MfUsgWel
-from flopy.modflow import ModflowBas, ModflowChd, ModflowDis
+from flopy.mfusg import (
+    MfUsg,
+    MfUsgBas,
+    MfUsgBcf,
+    MfUsgBct,
+    MfUsgCln,
+    MfUsgOc,
+    MfUsgSms,
+    MfUsgWel,
+)
+from flopy.modflow import ModflowChd, ModflowDis
 from flopy.plot import PlotCrossSection, PlotMapView
 from flopy.utils import HeadFile
 
@@ -98,7 +107,7 @@
 strt[1, :, 0] = 60.0
 strt[1, :, -1] = 60.0
 
-bas = ModflowBas(mf, ibound=ibound, strt=strt)
+bas = MfUsgBas(mf, ibound=ibound, strt=strt)
 bas.ibound.fmtin = "(25I3)"
 bas.strt.fmtin = "(10e12.4)"
 # -
diff --git a/.docs/Notebooks/mfusg_transport_Ex4_DualDomain.py b/.docs/Notebooks/mfusg_transport_tutorial04_DualDomain.py
similarity index 96%
rename from .docs/Notebooks/mfusg_transport_Ex4_DualDomain.py
rename to .docs/Notebooks/mfusg_transport_tutorial04_DualDomain.py
index 92ed5e979e..b0b5d264d7 100644
--- a/.docs/Notebooks/mfusg_transport_Ex4_DualDomain.py
+++ b/.docs/Notebooks/mfusg_transport_tutorial04_DualDomain.py
@@ -63,8 +63,17 @@
 import numpy as np
 
 import flopy
-from flopy.mfusg import MfUsg, MfUsgBcf, MfUsgBct, MfUsgDpt, MfUsgOc, MfUsgPcb, MfUsgSms
-from flopy.modflow import ModflowBas, ModflowChd, ModflowDis
+from flopy.mfusg import (
+    MfUsg,
+    MfUsgBas,
+    MfUsgBcf,
+    MfUsgBct,
+    MfUsgDpt,
+    MfUsgOc,
+    MfUsgPcb,
+    MfUsgSms,
+)
+from flopy.modflow import ModflowChd, ModflowDis
 from flopy.utils import HeadFile
 
 # -
@@ -121,7 +130,7 @@
 strt[:, :, 0] = 10.0
 strt[:, :, -1] = 9.0
 
-bas = ModflowBas(mf, ibound=ibound, strt=strt)
+bas = MfUsgBas(mf, ibound=ibound, strt=strt)
 # -
 # +
 tran = 1000.0
diff --git a/.docs/Notebooks/mfusg_transport_Ex5_Henry.py b/.docs/Notebooks/mfusg_transport_tutorial05_Henry.py
similarity index 98%
rename from .docs/Notebooks/mfusg_transport_Ex5_Henry.py
rename to .docs/Notebooks/mfusg_transport_tutorial05_Henry.py
index 69c5153f4b..c453992602 100644
--- a/.docs/Notebooks/mfusg_transport_Ex5_Henry.py
+++ b/.docs/Notebooks/mfusg_transport_tutorial05_Henry.py
@@ -33,6 +33,7 @@
 import flopy
 from flopy.mfusg import (
     MfUsg,
+    MfUsgBas,
     MfUsgBct,
     MfUsgDdf,
     MfUsgDisU,
@@ -42,7 +43,7 @@
     MfUsgSms,
     MfUsgWel,
 )
-from flopy.modflow import ModflowBas, ModflowChd, ModflowDis
+from flopy.modflow import ModflowChd, ModflowDis
 from flopy.plot import PlotCrossSection
 from flopy.utils import HeadUFile
 from flopy.utils.gridgen import Gridgen
@@ -52,7 +53,7 @@
 model_ws = "Ex5_Henry"
 mf = MfUsg(
     version="mfusg",
-    structured=True,
+    structured=False,
     model_ws=model_ws,
     modelname="Ex5_Henry",
     exe_name="mfusg_gsi",
@@ -91,7 +92,7 @@
 # mf.modelgrid = ugrid
 # -
 # +
-bas = ModflowBas(mf)
+bas = MfUsgBas(mf)
 # -
 # +
 ipakcb = 53
diff --git a/.docs/Notebooks/mfusg_transport_Ex6_Stallman_Heat_Transport.py b/.docs/Notebooks/mfusg_transport_tutorial06_StallmanHeatTransport.py
similarity index 88%
rename from .docs/Notebooks/mfusg_transport_Ex6_Stallman_Heat_Transport.py
rename to .docs/Notebooks/mfusg_transport_tutorial06_StallmanHeatTransport.py
index a958b93041..0c659dc834 100644
--- a/.docs/Notebooks/mfusg_transport_Ex6_Stallman_Heat_Transport.py
+++ b/.docs/Notebooks/mfusg_transport_tutorial06_StallmanHeatTransport.py
@@ -33,6 +33,7 @@
 import flopy
 from flopy.mfusg import (
     MfUsg,
+    MfUsgBas,
     MfUsgBct,
     MfUsgDdf,
     MfUsgDisU,
@@ -42,7 +43,7 @@
     MfUsgSms,
     MfUsgWel,
 )
-from flopy.modflow import ModflowBas, ModflowChd, ModflowDis
+from flopy.modflow import ModflowChd, ModflowDis
 from flopy.plot import PlotCrossSection
 from flopy.utils import HeadUFile
 from flopy.utils.gridgen import Gridgen
@@ -52,7 +53,7 @@
 model_ws = "Ex6a_Stallman_Heat"
 mf = MfUsg(
     version="mfusg",
-    structured=True,
+    structured=False,
     model_ws=model_ws,
     modelname="Ex6_Stallman",
     exe_name="mfusg_gsi",
@@ -103,7 +104,7 @@
 # -
 # +
 
-bas = ModflowBas(mf, strt=60.0)
+bas = MfUsgBas(mf, strt=60.0)
 # -
 
 # +
@@ -310,29 +311,29 @@
 success, buff = mf.run_model(silent=True)
 # -
 # +
-# vartype = [
-#     ("kstp", "<i4"),
-#     ("kper", "<i4"),
-#     ("pertim", "<f4"),
-#     ("totim", "<f4"),
-#     ("text", "S16"),
-#     ("nstrt", "<i4"),
-#     ("nndlay", "<i4"),
-#     ("ilay", "<i4"),
-#     ("data", "<f4"),
-# ]
-# hdr = np.fromfile(f"{mf.model_ws}/{mf.name}.con", vartype)
+vartype = [
+    ("kstp", "<i4"),
+    ("kper", "<i4"),
+    ("pertim", "<f4"),
+    ("totim", "<f4"),
+    ("text", "S16"),
+    ("nstrt", "<i4"),
+    ("nndlay", "<i4"),
+    ("ilay", "<i4"),
+    ("data", "<f4"),
+]
+hdr = np.fromfile(f"{mf.model_ws}/{mf.name}.con", vartype)
 
-# tmpr_data = hdr[np.char.strip(hdr['text'].astype(str)) == 'TMPR']
-# conc_data = hdr[np.char.strip(hdr['text'].astype(str)) == 'CONC']
-# simtemp = tmpr_data[(tmpr_data['kper']==601) & (tmpr_data['kstp']==6)]['data']
-# simconc = conc_data[(conc_data['kper']==601) & (conc_data['kstp']==6)]['data']
+tmpr_data = hdr[np.char.strip(hdr["text"].astype(str)) == "TMPR"]
+conc_data = hdr[np.char.strip(hdr["text"].astype(str)) == "CONC"]
+simtemp = tmpr_data[(tmpr_data["kper"] == 601) & (tmpr_data["kstp"] == 6)]["data"]
+simconc = conc_data[(conc_data["kper"] == 601) & (conc_data["kstp"] == 6)]["data"]
 
-tempobj = HeadUFile(f"{mf.model_ws}/{mf.name}.con", text="tmpr")
-simtemp = np.concatenate(tempobj.get_data())
-
-concobj = HeadUFile(f"{mf.model_ws}/{mf.name}.con", text="conc")
-simconc = concobj.get_data()
+# tempobj = HeadUFile(f"{mf.model_ws}/{mf.name}.con", precision="single", text=' TMPR           ')
+# simtemp = np.concatenate(tempobj.get_data())
+#
+# concobj = HeadUFile(f"{mf.model_ws}/{mf.name}.con", text="CONC")
+# simconc = concobj.get_data()
 # -
 # +
 depth = top - np.squeeze(dis.botm.array)
diff --git a/.docs/Notebooks/mfusg_transport_Ex7a_Discrete_Fracture.py b/.docs/Notebooks/mfusg_transport_tutorial07a_DiscreteFracture.py
similarity index 99%
rename from .docs/Notebooks/mfusg_transport_Ex7a_Discrete_Fracture.py
rename to .docs/Notebooks/mfusg_transport_tutorial07a_DiscreteFracture.py
index 2df6c9f3a6..c9c3bea966 100644
--- a/.docs/Notebooks/mfusg_transport_Ex7a_Discrete_Fracture.py
+++ b/.docs/Notebooks/mfusg_transport_tutorial07a_DiscreteFracture.py
@@ -35,6 +35,7 @@
 import flopy.utils.binaryfile as bf
 from flopy.mfusg import (
     MfUsg,
+    MfUsgBas,
     MfUsgBct,
     MfUsgDisU,
     MfUsgLpf,
@@ -44,7 +45,7 @@
     MfUsgSms,
     MfUsgWel,
 )
-from flopy.modflow import ModflowBas, ModflowChd, ModflowDis
+from flopy.modflow import ModflowChd, ModflowDis
 from flopy.plot import PlotCrossSection, PlotMapView
 from flopy.utils import HeadUFile
 from flopy.utils.gridgen import Gridgen
@@ -66,7 +67,7 @@
 
 mf = MfUsg(
     version="mfusg",
-    structured=True,
+    structured=False,
     model_ws=model_ws,
     modelname="Ex7_MD",
     exe_name="mfusg_gsi",
@@ -107,7 +108,7 @@
 disu = MfUsgDisU(mf, **gridprops, itmuni=5, lenuni=1, nper=nper, perlen=perlen)
 # -
 # +
-bas = ModflowBas(mf, strt=10.0)
+bas = MfUsgBas(mf, strt=10.0)
 # -
 # +
 
diff --git a/.docs/Notebooks/mfusg_transport_Ex7b_SandBox.py b/.docs/Notebooks/mfusg_transport_tutorial07b_SandBox.py
similarity index 99%
rename from .docs/Notebooks/mfusg_transport_Ex7b_SandBox.py
rename to .docs/Notebooks/mfusg_transport_tutorial07b_SandBox.py
index 470f73daf9..dcac8c3eca 100644
--- a/.docs/Notebooks/mfusg_transport_Ex7b_SandBox.py
+++ b/.docs/Notebooks/mfusg_transport_tutorial07b_SandBox.py
@@ -68,7 +68,7 @@
 model_ws = "Ex7_SandBox"
 mf = MfUsg(
     version="mfusg",
-    structured=True,
+    structured=False,
     model_ws=model_ws,
     modelname="Ex7_SandBox",
     exe_name="mfusg_gsi",
diff --git a/.docs/Notebooks/mfusg_transport_Ex7c_ChainDecay.py b/.docs/Notebooks/mfusg_transport_tutorial07c_ChainDecay.py
similarity index 98%
rename from .docs/Notebooks/mfusg_transport_Ex7c_ChainDecay.py
rename to .docs/Notebooks/mfusg_transport_tutorial07c_ChainDecay.py
index d8c4635f14..fe3b4fb4c9 100644
--- a/.docs/Notebooks/mfusg_transport_Ex7c_ChainDecay.py
+++ b/.docs/Notebooks/mfusg_transport_tutorial07c_ChainDecay.py
@@ -33,6 +33,7 @@
 import flopy.utils.binaryfile as bf
 from flopy.mfusg import (
     MfUsg,
+    MfUsgBas,
     MfUsgBct,
     MfUsgDisU,
     MfUsgLpf,
@@ -42,7 +43,7 @@
     MfUsgSms,
     MfUsgWel,
 )
-from flopy.modflow import ModflowBas, ModflowChd, ModflowDis
+from flopy.modflow import ModflowChd, ModflowDis
 from flopy.plot import PlotCrossSection, PlotMapView
 from flopy.utils import HeadUFile
 from flopy.utils.gridgen import Gridgen
@@ -62,7 +63,7 @@
 
 mf = MfUsg(
     version="mfusg",
-    structured=True,
+    structured=False,
     model_ws=model_ws,
     modelname="Ex7_PCE",
     exe_name="mfusg_gsi",
@@ -117,7 +118,7 @@
 # -
 
 # +
-bas = ModflowBas(mf, strt=20.0)
+bas = MfUsgBas(mf, strt=20.0)
 # -
 
 # +
diff --git a/.docs/Notebooks/mfusg_transport_Ex8_Lake.py b/.docs/Notebooks/mfusg_transport_tutorial08_Lake.py
similarity index 79%
rename from .docs/Notebooks/mfusg_transport_Ex8_Lake.py
rename to .docs/Notebooks/mfusg_transport_tutorial08_Lake.py
index eddbb170bb..262c475803 100644
--- a/.docs/Notebooks/mfusg_transport_Ex8_Lake.py
+++ b/.docs/Notebooks/mfusg_transport_tutorial08_Lake.py
@@ -33,6 +33,7 @@
 
 # +
 import os
+from pathlib import Path
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -41,6 +42,7 @@
 import flopy.utils.binaryfile as bf
 from flopy.mfusg import (
     MfUsg,
+    MfUsgBas,
     MfUsgBcf,
     MfUsgBct,
     MfUsgDisU,
@@ -51,7 +53,7 @@
     MfUsgRch,
     MfUsgSms,
 )
-from flopy.modflow import ModflowBas, ModflowChd, ModflowDis, ModflowFhb, ModflowGage
+from flopy.modflow import ModflowChd, ModflowDis, ModflowFhb, ModflowGage
 from flopy.plot import PlotCrossSection, PlotMapView
 from flopy.utils import HeadUFile
 from flopy.utils.gridgen import Gridgen
@@ -63,7 +65,7 @@
 
 mf = MfUsg(
     version="mfusg",
-    structured=True,
+    structured=False,
     model_ws=model_ws,
     modelname="Ex8_Lake",
     exe_name="mfusg",
@@ -157,7 +159,7 @@
         laybnd[7:10, 7:10] = 0
     laybnd = laybnd.reshape(-1)
     ibound.append(laybnd)
-bas = ModflowBas(mf, ibound=ibound, strt=115.0)
+bas = MfUsgBas(mf, ibound=ibound, strt=115.0)
 bas.ibound.fmtin = "(17I4)"
 # -
 
@@ -322,7 +324,7 @@
 
 # +
 gages = [[-1, -37, 3]]
-files = [f"{model_ws}.got"]
+files = [f"{mf.name}.got"]
 gage = ModflowGage(mf, numgage=1, gage_data=gages, files=files)
 
 # -
@@ -344,19 +346,221 @@
 # -
 
 # +
-from autotest.conftest import get_example_data_path
 
-txtfile = (
-    get_example_data_path() / "mfusg_transport" / "Ex8_Lake" / "Analytical_Soln.txt"
-)
-soln = np.genfromtxt(txtfile, skip_header=2)
+time = [
+    0.00e00,
+    1.60e01,
+    3.23e01,
+    4.90e01,
+    6.60e01,
+    8.33e01,
+    1.01e02,
+    1.19e02,
+    1.37e02,
+    1.56e02,
+    1.75e02,
+    1.95e02,
+    2.15e02,
+    2.35e02,
+    2.56e02,
+    2.77e02,
+    2.98e02,
+    3.20e02,
+    3.43e02,
+    3.66e02,
+    3.89e02,
+    4.13e02,
+    4.37e02,
+    4.62e02,
+    4.87e02,
+    5.13e02,
+    5.39e02,
+    5.66e02,
+    5.93e02,
+    6.21e02,
+    6.50e02,
+    6.79e02,
+    7.08e02,
+    7.38e02,
+    7.69e02,
+    8.01e02,
+    8.33e02,
+    8.65e02,
+    8.99e02,
+    9.33e02,
+    9.67e02,
+    1.00e03,
+    1.04e03,
+    1.08e03,
+    1.11e03,
+    1.15e03,
+    1.19e03,
+    1.23e03,
+    1.27e03,
+    1.31e03,
+    1.35e03,
+    1.40e03,
+    1.44e03,
+    1.49e03,
+    1.53e03,
+    1.58e03,
+    1.63e03,
+    1.67e03,
+    1.72e03,
+    1.77e03,
+    1.83e03,
+    1.88e03,
+    1.93e03,
+    1.99e03,
+    2.04e03,
+    2.10e03,
+    2.16e03,
+    2.22e03,
+    2.28e03,
+    2.34e03,
+    2.40e03,
+    2.47e03,
+    2.53e03,
+    2.60e03,
+    2.67e03,
+    2.74e03,
+    2.81e03,
+    2.88e03,
+    2.95e03,
+    3.03e03,
+    3.10e03,
+    3.18e03,
+    3.26e03,
+    3.34e03,
+    3.42e03,
+    3.51e03,
+    3.60e03,
+    3.68e03,
+    3.77e03,
+    3.86e03,
+    3.96e03,
+    4.05e03,
+    4.15e03,
+    4.25e03,
+    4.35e03,
+    4.45e03,
+    4.56e03,
+    4.67e03,
+    4.77e03,
+    4.89e03,
+    5.00e03,
+]
+manual_calculation = [
+    100.000,
+    96.6139,
+    93.3794,
+    90.1930,
+    87.0376,
+    83.9044,
+    80.7895,
+    77.6918,
+    74.6126,
+    71.5545,
+    68.5221,
+    65.5215,
+    62.5610,
+    59.6505,
+    56.8017,
+    54.0276,
+    51.3421,
+    48.7593,
+    46.2923,
+    43.9523,
+    41.7488,
+    39.6872,
+    37.7690,
+    35.9925,
+    34.3533,
+    32.8457,
+    31.4636,
+    30.2008,
+    29.0501,
+    28.0026,
+    27.0492,
+    26.1808,
+    25.3893,
+    24.6672,
+    24.0078,
+    23.4050,
+    22.8535,
+    22.3485,
+    21.8857,
+    21.4614,
+    21.0721,
+    20.7148,
+    20.3868,
+    20.0856,
+    19.8090,
+    19.5550,
+    19.3218,
+    19.1075,
+    18.9085,
+    18.7227,
+    18.5489,
+    18.3852,
+    18.2304,
+    18.0836,
+    17.9438,
+    17.8102,
+    17.6823,
+    17.5596,
+    17.4411,
+    17.3260,
+    17.2140,
+    17.1047,
+    16.9978,
+    16.8929,
+    16.7899,
+    16.6884,
+    16.5882,
+    16.4888,
+    16.3900,
+    16.2914,
+    16.1930,
+    16.0946,
+    15.9961,
+    15.8974,
+    15.7982,
+    15.6986,
+    15.5984,
+    15.4976,
+    15.3960,
+    15.2936,
+    15.1904,
+    15.0862,
+    14.9810,
+    14.8747,
+    14.7674,
+    14.6589,
+    14.5493,
+    14.4386,
+    14.3266,
+    14.2134,
+    14.0990,
+    13.9833,
+    13.8663,
+    13.7481,
+    13.6287,
+    13.5079,
+    13.3859,
+    13.2626,
+    13.1381,
+    13.0123,
+    12.8852,
+]
+
 # -
 
 # +
 fig = plt.figure(figsize=(8, 5), dpi=150)
 ax = fig.add_subplot(111)
 ax.plot(got[:, 0], got[:, 3])
-ax.plot(soln[:, 0], soln[:, 1], linestyle="--")
+ax.plot(time, manual_calculation, linestyle="--")
 
 ax2 = ax.twinx()
 ax2.plot(got[:, 0], got[:, 2], color="b")
diff --git a/.docs/Notebooks/mfusg_transport_Ex9_PFAS.py b/.docs/Notebooks/mfusg_transport_tutorial09_PFAS.py
similarity index 93%
rename from .docs/Notebooks/mfusg_transport_Ex9_PFAS.py
rename to .docs/Notebooks/mfusg_transport_tutorial09_PFAS.py
index 472e3dd2a0..67726238dc 100644
--- a/.docs/Notebooks/mfusg_transport_Ex9_PFAS.py
+++ b/.docs/Notebooks/mfusg_transport_tutorial09_PFAS.py
@@ -43,6 +43,7 @@
 import flopy.utils.binaryfile as bf
 from flopy.mfusg import (
     MfUsg,
+    MfUsgBas,
     MfUsgBct,
     MfUsgDisU,
     MfUsgLpf,
@@ -50,7 +51,7 @@
     MfUsgRch,
     MfUsgSms,
 )
-from flopy.modflow import ModflowBas, ModflowChd, ModflowDis
+from flopy.modflow import ModflowChd, ModflowDis
 from flopy.plot import PlotCrossSection, PlotMapView
 from flopy.utils import HeadUFile
 from flopy.utils.gridgen import Gridgen
@@ -61,7 +62,7 @@
 model_ws = "Ex9_PFAS"
 mf = MfUsg(
     version="mfusg",
-    structured=True,
+    structured=False,
     model_ws=model_ws,
     modelname="Ex9_PFAS",
     exe_name="mfusg_gsi",
@@ -92,26 +93,14 @@
     os.mkdir(gridgen_ws)
 g = Gridgen(ms.modelgrid, model_ws=gridgen_ws)
 g.build()
-# -
-# +
-disu = g.get_disu(mf, itmuni=3, lenuni=3, nper=1, perlen=2.025)
-disu.ivsd = -1
+gridprops = g.get_gridprops_disu5()
 anglex = g.get_anglex()
-disu.iac.fmtin = "(10I4)"
-disu.ja.fmtin = "(10I4)"
-disu.cl12.fmtin = "(10F6.2)"
-disu.fahl.fmtin = "(10F6.2)"
 # -
-# MODFLOW-USG does not have vertices, so we need to create
-# and unstructured grid and then assign it to the model. This
-# will allow plotting and other features to work properly.
 # +
-gridprops_ug = g.get_gridprops_unstructuredgrid()
-ugrid = flopy.discretization.UnstructuredGrid(**gridprops_ug)
-mf.modelgrid = ugrid
+disu = MfUsgDisU(mf, **gridprops)
 # -
 # +
-bas = ModflowBas(mf, ibound=1, strt=15.0, richards=True, unstructured=True)
+bas = MfUsgBas(mf, ibound=1, strt=15.0, richards=True)
 # -
 # +
 ipakcb = 50
@@ -224,10 +213,10 @@
 # -
 # +
 mf.write_input()
-success, buff = mf.run_model()
+success, buff = mf.run_model(silent=False)
 # -
 # +
-concobj = HeadUFile(f"{mf.model_ws}/{mf.name}.con", text="conc")
+concobj = HeadUFile(f"{mf.model_ws}/{mf.name}.con", text="CONC")
 simconc1 = concobj.get_ts((119))
 # -
 
diff --git a/.docs/tutorials.rst b/.docs/tutorials.rst
index 0173716121..8f3e6d318d 100644
--- a/.docs/tutorials.rst
+++ b/.docs/tutorials.rst
@@ -65,6 +65,23 @@ MODFLOW-NWT
    Notebooks/nwt_option_blocks_tutorial
 
 
+MODFLOW-USG
+-----------
+
+.. toctree::
+   :maxdepth: 2
+
+   Notebooks/mfusg_transport_tutorial01_1D
+   Notebooks/mfusg_transport_tutorial02_Radial
+   Notebooks/mfusg_transport_tutorial03_Conduit
+   Notebooks/mfusg_transport_tutorial04_DualDomain
+   Notebooks/mfusg_transport_tutorial05_Henry
+   Notebooks/mfusg_transport_tutorial06_StallmanHeatTransport
+   Notebooks/mfusg_transport_tutorial07a_DiscreteFracture
+   Notebooks/mfusg_transport_tutorial07b_SandBox
+   Notebooks/mfusg_transport_tutorial09_PFAS
+
+
 MT3DMS
 ------
 
diff --git a/autotest/test_example_notebooks.py b/autotest/test_example_notebooks.py
index 3855f8c19e..e7c673adf7 100644
--- a/autotest/test_example_notebooks.py
+++ b/autotest/test_example_notebooks.py
@@ -10,6 +10,8 @@
 
 EXCLUDE = [
     "mf6_lgr",
+    "mfusg_transport_tutorial07c_ChainDecay",
+    "mfusg_transport_tutorial08_Lake",
 ]
 # skip pyvista notebooks on windows/mac in CI due to persistent issues
 # first with offscreen rendering, then with finding pyvista even after
diff --git a/examples/data/mfusg_transport/Ex8_Lake/Analytical_Soln.txt b/examples/data/mfusg_transport/Ex8_Lake/Analytical_Soln.txt
deleted file mode 100644
index 25708afdf7..0000000000
--- a/examples/data/mfusg_transport/Ex8_Lake/Analytical_Soln.txt
+++ /dev/null
@@ -1,103 +0,0 @@
-
-Time        Manual_Calculation
-0.00E+00    100.000
-1.60E+01    96.6139
-3.23E+01    93.3794
-4.90E+01    90.1930
-6.60E+01    87.0376
-8.33E+01    83.9044
-1.01E+02    80.7895
-1.19E+02    77.6918
-1.37E+02    74.6126
-1.56E+02    71.5545
-1.75E+02    68.5221
-1.95E+02    65.5215
-2.15E+02    62.5610
-2.35E+02    59.6505
-2.56E+02    56.8017
-2.77E+02    54.0276
-2.98E+02    51.3421
-3.20E+02    48.7593
-3.43E+02    46.2923
-3.66E+02    43.9523
-3.89E+02    41.7488
-4.13E+02    39.6872
-4.37E+02    37.7690
-4.62E+02    35.9925
-4.87E+02    34.3533
-5.13E+02    32.8457
-5.39E+02    31.4636
-5.66E+02    30.2008
-5.93E+02    29.0501
-6.21E+02    28.0026
-6.50E+02    27.0492
-6.79E+02    26.1808
-7.08E+02    25.3893
-7.38E+02    24.6672
-7.69E+02    24.0078
-8.01E+02    23.4050
-8.33E+02    22.8535
-8.65E+02    22.3485
-8.99E+02    21.8857
-9.33E+02    21.4614
-9.67E+02    21.0721
-1.00E+03    20.7148
-1.04E+03    20.3868
-1.08E+03    20.0856
-1.11E+03    19.8090
-1.15E+03    19.5550
-1.19E+03    19.3218
-1.23E+03    19.1075
-1.27E+03    18.9085
-1.31E+03    18.7227
-1.35E+03    18.5489
-1.40E+03    18.3852
-1.44E+03    18.2304
-1.49E+03    18.0836
-1.53E+03    17.9438
-1.58E+03    17.8102
-1.63E+03    17.6823
-1.67E+03    17.5596
-1.72E+03    17.4411
-1.77E+03    17.3260
-1.83E+03    17.2140
-1.88E+03    17.1047
-1.93E+03    16.9978
-1.99E+03    16.8929
-2.04E+03    16.7899
-2.10E+03    16.6884
-2.16E+03    16.5882
-2.22E+03    16.4888
-2.28E+03    16.3900
-2.34E+03    16.2914
-2.40E+03    16.1930
-2.47E+03    16.0946
-2.53E+03    15.9961
-2.60E+03    15.8974
-2.67E+03    15.7982
-2.74E+03    15.6986
-2.81E+03    15.5984
-2.88E+03    15.4976
-2.95E+03    15.3960
-3.03E+03    15.2936
-3.10E+03    15.1904
-3.18E+03    15.0862
-3.26E+03    14.9810
-3.34E+03    14.8747
-3.42E+03    14.7674
-3.51E+03    14.6589
-3.60E+03    14.5493
-3.68E+03    14.4386
-3.77E+03    14.3266
-3.86E+03    14.2134
-3.96E+03    14.0990
-4.05E+03    13.9833
-4.15E+03    13.8663
-4.25E+03    13.7481
-4.35E+03    13.6287
-4.45E+03    13.5079
-4.56E+03    13.3859
-4.67E+03    13.2626
-4.77E+03    13.1381
-4.89E+03    13.0123
-5.00E+03    12.8852
\ No newline at end of file
diff --git a/flopy/mfusg/mfusgbas.py b/flopy/mfusg/mfusgbas.py
index b2cf958b36..544eefc46f 100644
--- a/flopy/mfusg/mfusgbas.py
+++ b/flopy/mfusg/mfusgbas.py
@@ -38,6 +38,9 @@ class MfUsgBas(Package):
     ichflg : bool, optional
         Flag indicating that flows between constant head cells should be
         calculated (the default is False).
+    richards: bool, optional
+        Flag indicating that Richards equation is being solved (the default
+        is False). This is equivalent to using laytyp=5 in the LPF package.
     stoper : float
         percent discrepancy that is compared to the budget percent discrepancy
         continue when the solver convergence criteria are not met.  Execution
diff --git a/flopy/mfusg/mfusgbct.py b/flopy/mfusg/mfusgbct.py
index 93de791786..b2fd3a11a0 100644
--- a/flopy/mfusg/mfusgbct.py
+++ b/flopy/mfusg/mfusgbct.py
@@ -459,7 +459,7 @@ def write_file(self, f=None):
         if f is None:
             f_obj = open(self.fn_path, "w")
 
-        # Item 1: ITRNSP ipakcb MCOMP ICBNDFLG ITVD IADSORB ICT CINACT CICLOSE
+        # Item 1a: ITRNSP ipakcb MCOMP ICBNDFLG ITVD IADSORB ICT CINACT CICLOSE
         # IDISP IXDISP DIFFNC IZOD IFOD IFMBC IHEAT IMCOMP IDISPCLN NSEQITR
         f_obj.write(f"{self.heading}\n")
 
@@ -488,6 +488,7 @@ def write_file(self, f=None):
 
         f_obj.write(self.options + "\n")
 
+        # Item 1b.
         if self.ifmbc:
             f_obj.write(
                 f" {self.unit_number[1]:9d} {self.unit_number[2]:9d}"
diff --git a/flopy/mfusg/mfusglpf.py b/flopy/mfusg/mfusglpf.py
index a965a41389..b8f5e69166 100644
--- a/flopy/mfusg/mfusglpf.py
+++ b/flopy/mfusg/mfusglpf.py
@@ -44,11 +44,13 @@ class MfUsgLpf(ModflowLpf):
         (default is -1.e30).
     laytyp : int or array of ints (nlay)
         Layer type, contains a flag for each layer that specifies the layer
-        type.
-        0 confined
-        >0 convertible
-        <0 convertible unless the THICKSTRT option is in effect.
-        (default is 0).
+        type (default is 0).
+        0 is confined
+        >0 is convertible
+        <0 is convertible unless the THICKSTRT option is in effect.
+        =4 is convertible, with transmissivity computed using upstream
+        water-table depth.
+        =5 is richards equation for variably saturated flow.
     layavg : int or array of ints (nlay)
         Layer average
         0 is harmonic mean
@@ -123,7 +125,7 @@ class MfUsgLpf(ModflowLpf):
         When STORAGECOEFFICIENT is used, Ss is confined storage coefficient.
         (default is 1.e-5).
     sy : float or array of floats (nlay, nrow, ncol)
-        is specific yield.
+        is specific yield or the saturated water content if richards is used.
         (default is 0.15).
     vkcb : float or array of floats (nlay, nrow, ncol)
         is the vertical hydraulic conductivity of a Quasi-three-dimensional
@@ -160,7 +162,20 @@ class MfUsgLpf(ModflowLpf):
          5-8 of USGS Techniques and Methods Report 6-A16 and the vertical
          conductance correction described on p. 5-18 of that report.
          (default is False).
-    extension : string
+    bubblept : boolean
+        turns on the use of the bubbling point pressure in the unsaturated zone
+        calculation such that gas saturations are represented.
+    alpha : float or array of floats (nlay, nrow, ncol)
+        is the van Genuchten shape parameter alpha
+    beta : float or array of floats (nlay, nrow, ncol)
+        is the van Genuchten shape parameter beta (also called n sometimes)
+    sr : float or array of floats (nlay, nrow, ncol)
+        is the residual saturation (equal to theta_r/theta_s)
+    brook : float or array of floats (nlay, nrow, ncol)
+        is the Brooks-Corey shape parameter for the hydraulic conductivity function
+    bp : float or array of floats (nlay, nrow, ncol)
+        is the bubbling point pressure / air entry head. Requires bubblept to be True.
+        extension : string
         Filename extension (default is 'lpf')
     unitnumber : int
         File unit number (default is None).
@@ -233,11 +248,11 @@ def __init__(
         novfc=False,
         bubblept=False,
         fullydry=False,
-        alpha=0,
-        beta=0,
-        sr=0,
-        brook=0,
-        bp=0,
+        alpha=1.0,
+        beta=7.0,
+        sr=0.05,
+        brook=6.0,
+        bp=0.0,
         extension="lpf",
         unitnumber=None,
         filenames=None,
@@ -320,8 +335,17 @@ def __init__(
                 model, (njag,), np.float32, ksat, "ksat", locat=self.unit_number[0]
             )
 
-        bas = model.get_package("BAS6")
-        self.richards = bas.richards
+        if self.laytyp == 5:
+            self.richards = True
+            bas = model.get_package("BAS6")
+            if not hasattr(bas, "richards") or not bas.richards:
+                raise ValueError(
+                    "The MfUsgBas package must have richards=True"
+                    "when using laytyp=5 in the LPF package."
+                )
+        else:
+            self.richards = False
+
         self.bubblept = bubblept
         self.fullydry = fullydry
 
diff --git a/flopy/utils/gridgen.py b/flopy/utils/gridgen.py
index 5401daded5..efe66c82c8 100644
--- a/flopy/utils/gridgen.py
+++ b/flopy/utils/gridgen.py
@@ -1225,6 +1225,34 @@ def get_angldegx(self, fldr=None):
         angldegx = np.where(fldr == -2, 270, angldegx)
         return angldegx
 
+    def get_anglex(self, fldr=None):
+        """
+        Get the anglex array
+
+        Parameters
+        ----------
+        fldr : ndarray
+            Flow direction indicator array.  If None, then it is read from
+            gridgen output.
+
+        Returns
+        -------
+        anglex : ndarray
+            A 1D vector indicating the angle (in radians) between the x
+            axis and an outward normal to the face.
+
+        """
+
+        if fldr is None:
+            fldr = self.get_fldr()
+        anglex = np.zeros(fldr.shape, dtype=float)
+        anglex = np.where(fldr == 0, 0.0, anglex)
+        anglex = np.where(abs(fldr) == 3, 0.0, anglex)
+        anglex = np.where(fldr == -2, 1.570796, anglex)
+        anglex = np.where(fldr == -1, 3.141593, anglex)
+        anglex = np.where(fldr == 2, 4.712389, anglex)
+        return anglex
+
     def get_verts_iverts(self, ncells, verbose=False):
         """
         Return a 2d array of x and y vertices and a list of size ncells that