diff --git a/.build_rtd_docs/conf.py b/.build_rtd_docs/conf.py
index a29ea3e0ca1..1d4829f789e 100644
--- a/.build_rtd_docs/conf.py
+++ b/.build_rtd_docs/conf.py
@@ -1,3 +1,4 @@
+# ruff: noqa: E402
 # Configuration file for the Sphinx documentation builder.
 #
 # This file only contains a selection of the most common options. For a full
@@ -14,7 +15,7 @@
 # documentation root, use os.path.abspath to make it absolute, like shown here.
 #
 import sys
-from subprocess import Popen, PIPE
+from subprocess import PIPE, Popen
 
 sys.path.insert(0, os.path.abspath(os.path.join("..", "doc")))
 sys.path.insert(0, os.path.abspath(os.path.join("..", "distribution")))
@@ -37,6 +38,7 @@
 # -- Update the modflow 6 version -------------------------------------------
 print("Update the modflow6 version")
 from update_version import update_version
+
 update_version()
 
 # -- import version from doc/version.py -------------------------------------
diff --git a/.github/common/check_format.py b/.github/common/check_format.py
index f06a71be394..733909898f7 100644
--- a/.github/common/check_format.py
+++ b/.github/common/check_format.py
@@ -3,10 +3,10 @@
 import sys
 import timeit
 from itertools import repeat
-from queue import Empty
+from multiprocessing import Manager, Pool, cpu_count
 from pathlib import Path
+from queue import Empty
 from subprocess import run
-from multiprocessing import cpu_count, Pool, Manager
 
 PROJ_ROOT = Path(__file__).parents[2]
 
@@ -20,13 +20,11 @@
     PROJ_ROOT / "src" / "Utilities" / "Libraries" / "sparskit2",
     PROJ_ROOT / "utils" / "mf5to6",
     PROJ_ROOT / "utils" / "zonebudget" / "msvs",
-    PROJ_ROOT / "msvs"
+    PROJ_ROOT / "msvs",
 ]
 
 # exclude these files from checks
-excludefiles = [
-    PROJ_ROOT / "src" / "Idm" / "gwf-stoidm.f90"
-]
+excludefiles = [PROJ_ROOT / "src" / "Idm" / "gwf-stoidm.f90"]
 
 # commands
 fprettify = "fprettify -c .fprettify.yaml"
diff --git a/.github/common/check_vfproj.py b/.github/common/check_vfproj.py
index bd570c55fd7..6c1e8eade7c 100644
--- a/.github/common/check_vfproj.py
+++ b/.github/common/check_vfproj.py
@@ -4,7 +4,6 @@
 from pathlib import Path
 from pprint import pformat
 
-
 PROJ_ROOT = Path(__file__).parents[2]
 
 
diff --git a/.hpc/sstat_poll.py b/.hpc/sstat_poll.py
index 12e60f101a0..fa18c9e0df0 100644
--- a/.hpc/sstat_poll.py
+++ b/.hpc/sstat_poll.py
@@ -141,9 +141,7 @@ def _run_command(
         )
         is None
     ):
-        raise ValueError(
-            f"SLURM command '{slurm_args.command}' does not exist"
-        )
+        raise ValueError(f"SLURM command '{slurm_args.command}' does not exist")
 
     end_tag = f"sstat:,error:,no,steps,running,for,job,{slurm_args.jobid}"
     # open file
diff --git a/autotest/test_gwf_multimvr.py b/autotest/test_gwf_multimvr.py
index 8de11be4820..70f06827f53 100644
--- a/autotest/test_gwf_multimvr.py
+++ b/autotest/test_gwf_multimvr.py
@@ -1000,9 +1000,10 @@ def check_output(idx, test):
     # cur_ws, gwfparent = ex[idx], gwf_names[idx]
     cur_ws = test.workspace
     gwfparent = "gwf_" + cases[idx] + "_p"
-    with open(os.path.join(cur_ws, gwfparent + ".lst"), "r") as gwf_lst, open(
-        os.path.join(cur_ws, "mfsim.lst"), "r"
-    ) as sim_lst:
+    with (
+        open(os.path.join(cur_ws, gwfparent + ".lst"), "r") as gwf_lst,
+        open(os.path.join(cur_ws, "mfsim.lst"), "r") as sim_lst,
+    ):
         gwf_lst_lines = gwf_lst.readlines()
         sim_lst_lines = sim_lst.readlines()
 
diff --git a/doc/SuppTechInfo/python/DRN-Discharge-Scaling.ipynb b/doc/SuppTechInfo/python/DRN-Discharge-Scaling.ipynb
index ae101e3636a..6d158dec356 100644
--- a/doc/SuppTechInfo/python/DRN-Discharge-Scaling.ipynb
+++ b/doc/SuppTechInfo/python/DRN-Discharge-Scaling.ipynb
@@ -7,9 +7,10 @@
    "outputs": [],
    "source": [
     "import os\n",
-    "import numpy as np\n",
+    "\n",
+    "import matplotlib.gridspec as gridspec\n",
     "import matplotlib.pyplot as plt\n",
-    "import matplotlib.gridspec as gridspec"
+    "import numpy as np"
    ]
   },
   {
@@ -20,6 +21,7 @@
    "source": [
     "try:\n",
     "    import spnspecs\n",
+    "\n",
     "    spnspecs.set_graph_specifications()\n",
     "except:\n",
     "    spnspecs = None"
@@ -31,7 +33,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "figpth = '../Figures'\n",
+    "figpth = \"../Figures\"\n",
     "width = 6.8\n",
     "dpi = 300"
    ]
@@ -44,9 +46,9 @@
    "source": [
     "def no_smoothing(sat):\n",
     "    f = sat.copy()\n",
-    "    f[sat < 0] = 0.\n",
-    "    f[sat >= 0] = 1.\n",
-    "    return f    "
+    "    f[sat < 0] = 0.0\n",
+    "    f[sat >= 0] = 1.0\n",
+    "    return f"
    ]
   },
   {
@@ -57,9 +59,9 @@
    "source": [
     "def linear_smoothing(sat):\n",
     "    f = sat.copy()\n",
-    "    f[sat < 0] = 0.\n",
-    "    f[sat > 1] = 1.\n",
-    "    return f    "
+    "    f[sat < 0] = 0.0\n",
+    "    f[sat > 1] = 1.0\n",
+    "    return f"
    ]
   },
   {
@@ -82,7 +84,7 @@
    "outputs": [],
    "source": [
     "def cubic_smoothing_dervsat(sat, c1=-1, c2=2):\n",
-    "    dfdx = 3 * c1 * sat**2 + 2. * c2 * sat\n",
+    "    dfdx = 3 * c1 * sat**2 + 2.0 * c2 * sat\n",
     "    dfdx[sat < 0] = 0\n",
     "    dfdx[sat > 1] = 0\n",
     "    return dfdx"
@@ -96,8 +98,8 @@
    "source": [
     "def cubic_smoothingalt(x, xrange, c1=-1, c2=2):\n",
     "    sat = x / xrange\n",
-    "    cof1 = c1 / (xrange)**3\n",
-    "    cof2 = c2 / (xrange)**2\n",
+    "    cof1 = c1 / (xrange) ** 3\n",
+    "    cof2 = c2 / (xrange) ** 2\n",
     "    f = cof1 * x**3 + cof2 * x**2\n",
     "    f[sat < 0] = 0\n",
     "    f[sat > 1] = 1\n",
@@ -112,8 +114,8 @@
    "source": [
     "def cubic_smoothing_dervh(x, xrange, c1=-1, c2=2):\n",
     "    sat = x / xrange\n",
-    "    cof1 = 3 * c1 / (xrange)**3\n",
-    "    cof2 = 2 * c2 / (xrange)**2\n",
+    "    cof1 = 3 * c1 / (xrange) ** 3\n",
+    "    cof2 = 2 * c2 / (xrange) ** 2\n",
     "    dfdx = cof1 * x**2 + cof2 * x\n",
     "    dfdx[sat < 0] = 0\n",
     "    dfdx[sat > 1] = 0\n",
@@ -129,9 +131,9 @@
     "def fd_derv(x, v):\n",
     "    isize = v.shape[0]\n",
     "    derv = np.ones((isize), dtype=v.dtype) * np.nan\n",
-    "    for i in range(1, isize-2):\n",
-    "        dv = v[i+1] - v[i-1]\n",
-    "        dx = x[i+1] - x[i-1]\n",
+    "    for i in range(1, isize - 2):\n",
+    "        dv = v[i + 1] - v[i - 1]\n",
+    "        dx = x[i + 1] - x[i - 1]\n",
     "        derv[i] = dv / dx\n",
     "    return derv"
    ]
@@ -147,8 +149,8 @@
     "s0, s1 = -0.1, 1.1\n",
     "s = np.arange(s0, s1, step=step, dtype=dtype)\n",
     "sp = s.copy()\n",
-    "sp[s < 0] = 0.\n",
-    "sp[s > 1] = 1."
+    "sp[s < 0] = 0.0\n",
+    "sp[s > 1] = 1.0"
    ]
   },
   {
@@ -166,7 +168,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(s, sp*sp)"
+    "plt.plot(s, sp * sp)"
    ]
   },
   {
@@ -176,7 +178,7 @@
    "outputs": [],
    "source": [
     "f0 = cubic_smoothing(s)\n",
-    "f1 = cubic_smoothing(s, c1=-2., c2=3.)"
+    "f1 = cubic_smoothing(s, c1=-2.0, c2=3.0)"
    ]
   },
   {
@@ -185,11 +187,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "dseep = 1.\n",
-    "z0 = 0.\n",
+    "dseep = 1.0\n",
+    "z0 = 0.0\n",
     "xdiff = dseep * s\n",
     "h = xdiff + z0\n",
-    "dsdh = 1. / dseep"
+    "dsdh = 1.0 / dseep"
    ]
   },
   {
@@ -198,13 +200,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(s, sp, lw=1.5, color='black', label='Linear')\n",
-    "plt.plot(s, f0, lw=0.75, color='blue', label='GWSEEP')\n",
-    "plt.plot(s, f1, lw=0.75, color='red', label='WELL')\n",
+    "plt.plot(s, sp, lw=1.5, color=\"black\", label=\"Linear\")\n",
+    "plt.plot(s, f0, lw=0.75, color=\"blue\", label=\"GWSEEP\")\n",
+    "plt.plot(s, f1, lw=0.75, color=\"red\", label=\"WELL\")\n",
     "if spnspecs is not None:\n",
     "    spnspecs.graph_legend(plt.gca())\n",
-    "#plt.xlim(0.9, 1.1)\n",
-    "#plt.ylim(0.9, 1.1)"
+    "# plt.xlim(0.9, 1.1)\n",
+    "# plt.ylim(0.9, 1.1)"
    ]
   },
   {
@@ -214,7 +216,7 @@
    "outputs": [],
    "source": [
     "f0p = cubic_smoothingalt(xdiff, dseep)\n",
-    "f1p = cubic_smoothingalt(xdiff, dseep, c1=-2., c2=3.)"
+    "f1p = cubic_smoothingalt(xdiff, dseep, c1=-2.0, c2=3.0)"
    ]
   },
   {
@@ -223,10 +225,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(s, f0, lw=0.75, color='blue', label='GWSEEP')\n",
-    "plt.plot(s, f0p, lw=0., color='blue', marker='o', mfc='none', ms=4, markevery=50)\n",
-    "plt.plot(s, f1, lw=0.75, color='red', label='WELL')\n",
-    "plt.plot(s, f1p, lw=0., color='red', marker='o', mfc='none', ms=4, markevery=50)"
+    "plt.plot(s, f0, lw=0.75, color=\"blue\", label=\"GWSEEP\")\n",
+    "plt.plot(s, f0p, lw=0.0, color=\"blue\", marker=\"o\", mfc=\"none\", ms=4, markevery=50)\n",
+    "plt.plot(s, f1, lw=0.75, color=\"red\", label=\"WELL\")\n",
+    "plt.plot(s, f1p, lw=0.0, color=\"red\", marker=\"o\", mfc=\"none\", ms=4, markevery=50)"
    ]
   },
   {
@@ -235,8 +237,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(s, f0-f0p, lw=0.75, color='blue', label='GWSEEP')\n",
-    "plt.plot(s, f1-f1p, lw=0.75, color='red', label='WELL')"
+    "plt.plot(s, f0 - f0p, lw=0.75, color=\"blue\", label=\"GWSEEP\")\n",
+    "plt.plot(s, f1 - f1p, lw=0.75, color=\"red\", label=\"WELL\")"
    ]
   },
   {
@@ -246,8 +248,8 @@
    "outputs": [],
    "source": [
     "spderv = np.ones(s.shape, dtype=np.float)\n",
-    "spderv[s < 0] = 0.\n",
-    "spderv[s > 1] = 0."
+    "spderv[s < 0] = 0.0\n",
+    "spderv[s > 1] = 0.0"
    ]
   },
   {
@@ -257,7 +259,7 @@
    "outputs": [],
    "source": [
     "f0derv = cubic_smoothing_dervsat(s)\n",
-    "f1derv = cubic_smoothing_dervsat(s, c1=-2., c2=3.)"
+    "f1derv = cubic_smoothing_dervsat(s, c1=-2.0, c2=3.0)"
    ]
   },
   {
@@ -267,7 +269,7 @@
    "outputs": [],
    "source": [
     "f0derv2 = cubic_smoothing_dervh(xdiff, dseep)\n",
-    "f1derv2 = cubic_smoothing_dervh(xdiff, dseep, c1=-2., c2=3.)"
+    "f1derv2 = cubic_smoothing_dervh(xdiff, dseep, c1=-2.0, c2=3.0)"
    ]
   },
   {
@@ -276,10 +278,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(s, f0derv, lw=4, color='blue', label='GWSEEP', alpha=0.5)\n",
-    "plt.plot(s, f0derv2/dsdh, lw=0.75, color='black', ls='-.') #marker='o', mfc='none', ms=4, markevery=50)\n",
-    "plt.plot(s, f1derv, lw=4, color='red', label='WELL', alpha=0.5)\n",
-    "plt.plot(s, f1derv2/dsdh, lw=0.75, color='black', ls='-.') #marker='o', mfc='none', ms=4, markevery=50)"
+    "plt.plot(s, f0derv, lw=4, color=\"blue\", label=\"GWSEEP\", alpha=0.5)\n",
+    "plt.plot(\n",
+    "    s, f0derv2 / dsdh, lw=0.75, color=\"black\", ls=\"-.\"\n",
+    ")  # marker='o', mfc='none', ms=4, markevery=50)\n",
+    "plt.plot(s, f1derv, lw=4, color=\"red\", label=\"WELL\", alpha=0.5)\n",
+    "plt.plot(\n",
+    "    s, f1derv2 / dsdh, lw=0.75, color=\"black\", ls=\"-.\"\n",
+    ")  # marker='o', mfc='none', ms=4, markevery=50)"
    ]
   },
   {
@@ -288,9 +294,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(s, fd_derv(h, f0), lw=0.75, color='blue', label='GWSEEP')\n",
-    "plt.plot(s, fd_derv(h, f1), lw=0.75, color='red', label='WELL')\n",
-    "plt.ylabel(r'$\\frac{\\partial F}{\\partial h}$');"
+    "plt.plot(s, fd_derv(h, f0), lw=0.75, color=\"blue\", label=\"GWSEEP\")\n",
+    "plt.plot(s, fd_derv(h, f1), lw=0.75, color=\"red\", label=\"WELL\")\n",
+    "plt.ylabel(r\"$\\frac{\\partial F}{\\partial h}$\");"
    ]
   },
   {
@@ -299,9 +305,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(s, fd_derv(h, f0p), lw=0.75, color='blue', label='GWSEEP')\n",
-    "plt.plot(s, fd_derv(h, f1p), lw=0.75, color='red', label='WELL')\n",
-    "plt.ylabel(r'$\\frac{\\partial F}{\\partial h}$');"
+    "plt.plot(s, fd_derv(h, f0p), lw=0.75, color=\"blue\", label=\"GWSEEP\")\n",
+    "plt.plot(s, fd_derv(h, f1p), lw=0.75, color=\"red\", label=\"WELL\")\n",
+    "plt.ylabel(r\"$\\frac{\\partial F}{\\partial h}$\");"
    ]
   },
   {
@@ -324,46 +330,51 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig, axes = plt.subplots(nrows=1, ncols=2, tight_layout=True, figsize=(width, (1.4/s1) * width/2))\n",
+    "fig, axes = plt.subplots(\n",
+    "    nrows=1, ncols=2, tight_layout=True, figsize=(width, (1.4 / s1) * width / 2)\n",
+    ")\n",
     "\n",
-    "letters = ['A', 'B']\n",
+    "letters = [\"A\", \"B\"]\n",
     "for idx, ax in enumerate(axes):\n",
     "    ax.set_xlim(s0, s1)\n",
     "    ax.set_ylim(s0, 1.4)\n",
     "    ax.set_xticks([0, 0.25, 0.5, 0.75, 1])\n",
-    "    ax.set_xticklabels(['0', '0.25', '0.50', '0.75', '1.00'])\n",
+    "    ax.set_xticklabels([\"0\", \"0.25\", \"0.50\", \"0.75\", \"1.00\"])\n",
     "    ax.set_yticks([0, 0.25, 0.5, 0.75, 1, 1.25])\n",
-    "    ax.set_yticklabels(['0', '0.25', '0.50', '0.75', '1.00', '1.25'])\n",
+    "    ax.set_yticklabels([\"0\", \"0.25\", \"0.50\", \"0.75\", \"1.00\", \"1.25\"])\n",
     "    if spnspecs is not None:\n",
     "        spnspecs.remove_edge_ticks(ax)\n",
     "        spnspecs.heading(ax, letter=letters[idx])\n",
-    "    \n",
+    "\n",
     "ax = axes[0]\n",
-    "ax.axhline(0, lw=0.5, ls='-.', color='black')\n",
-    "ax.axhline(1, lw=0.5, ls='-.', color='black')\n",
-    "ax.axvline(0, lw=0.5, ls='-.', color='black')\n",
-    "ax.axvline(1, lw=0.5, ls='-.', color='black')\n",
-    "ax.plot(s, sp, lw=3.5, color='blue', label=r'$F_{DRN}$')\n",
-    "ax.plot(s, f0, lw=1.75, color='red', label=r'$F_{DRN}^*$')\n",
-    "ax.set_xlabel(r'$\\frac{h - ZDRN}{DDRN}$, unitless')\n",
-    "ax.set_ylabel('Discharge scale factor, unitless')\n",
+    "ax.axhline(0, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.axhline(1, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.axvline(0, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.axvline(1, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.plot(s, sp, lw=3.5, color=\"blue\", label=r\"$F_{DRN}$\")\n",
+    "ax.plot(s, f0, lw=1.75, color=\"red\", label=r\"$F_{DRN}^*$\")\n",
+    "ax.set_xlabel(r\"$\\frac{h - ZDRN}{DDRN}$, unitless\")\n",
+    "ax.set_ylabel(\"Discharge scale factor, unitless\")\n",
     "\n",
     "ax = axes[1]\n",
-    "ax.axhline(0, lw=0.5, ls='-.', color='black')\n",
-    "ax.axhline(1, lw=0.5, ls='-.', color='black')\n",
-    "ax.axvline(0, lw=0.5, ls='-.', color='black')\n",
-    "ax.axvline(1, lw=0.5, ls='-.', color='black')\n",
-    "ax.plot([0, 1], [1, 1], lw=3.5, color='blue', label=r'$F_{DRN}$ (linear)')\n",
-    "ax.plot([s0, 0], [0, 0], lw=3.5, color='blue', label=None)\n",
-    "ax.plot([1, s1], [0, 0], lw=3.5, color='blue', label=None)\n",
-    "ax.plot(s[s<=1], f0derv[s<=1], lw=1.75, color='red', label=r'$F_{DRN}^*$ (cubic)')\n",
-    "ax.plot([1, s1], [0, 0], lw=1.75, color='red', label=None)\n",
+    "ax.axhline(0, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.axhline(1, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.axvline(0, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.axvline(1, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.plot([0, 1], [1, 1], lw=3.5, color=\"blue\", label=r\"$F_{DRN}$ (linear)\")\n",
+    "ax.plot([s0, 0], [0, 0], lw=3.5, color=\"blue\", label=None)\n",
+    "ax.plot([1, s1], [0, 0], lw=3.5, color=\"blue\", label=None)\n",
+    "ax.plot(s[s <= 1], f0derv[s <= 1], lw=1.75, color=\"red\", label=r\"$F_{DRN}^*$ (cubic)\")\n",
+    "ax.plot([1, s1], [0, 0], lw=1.75, color=\"red\", label=None)\n",
     "if spnspecs is not None:\n",
-    "    spnspecs.graph_legend(ax, loc='lower right', bbox_to_anchor=(0.9,0.05))\n",
-    "ax.set_xlabel(r'$\\frac{h - ZDRN}{DDRN}$, unitless')\n",
-    "ax.set_ylabel('Discharge scale factor derivative with respect to ' + r'$\\frac{h - ZDRN}{DDRN}$, unitless')\n",
+    "    spnspecs.graph_legend(ax, loc=\"lower right\", bbox_to_anchor=(0.9, 0.05))\n",
+    "ax.set_xlabel(r\"$\\frac{h - ZDRN}{DDRN}$, unitless\")\n",
+    "ax.set_ylabel(\n",
+    "    \"Discharge scale factor derivative with respect to \"\n",
+    "    + r\"$\\frac{h - ZDRN}{DDRN}$, unitless\"\n",
+    ")\n",
     "\n",
-    "fpth = os.path.join(figpth, 'DischargeScaleFactor.pdf')\n",
+    "fpth = os.path.join(figpth, \"DischargeScaleFactor.pdf\")\n",
     "fig.savefig(fpth, dpi=dpi)"
    ]
   },
@@ -381,7 +392,7 @@
     "for idx in range(1, dc.shape[0]):\n",
     "    v = dc[idx]\n",
     "    if v > d:\n",
-    "        slope *= -1.\n",
+    "        slope *= -1.0\n",
     "        d += dadd\n",
     "    z[idx] = z[idx - 1] + slope * (v - dc[idx - 1])\n",
     "z -= z.max() / 2\n",
@@ -395,7 +406,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "t = np.arange(s.shape[0]) / s.shape[0] \n",
+    "t = np.arange(s.shape[0]) / s.shape[0]\n",
     "xdiffq = xdiff.copy()\n",
     "xdiffq[xdiff < 0] = 0\n",
     "q0 = no_smoothing(s) * xdiffq\n",
@@ -409,8 +420,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "arrowprops=dict(arrowstyle='->',\n",
-    "                connectionstyle=\"arc3\")"
+    "arrowprops = dict(arrowstyle=\"->\", connectionstyle=\"arc3\")"
    ]
   },
   {
@@ -419,98 +429,131 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig = plt.figure(tight_layout=True, figsize=(width, (2/3)*width))\n",
+    "fig = plt.figure(tight_layout=True, figsize=(width, (2 / 3) * width))\n",
     "gs = gridspec.GridSpec(2, 3)\n",
-    "axes = [fig.add_subplot(gs[0, :]), fig.add_subplot(gs[1, 0]),\n",
-    "        fig.add_subplot(gs[1, 1]), fig.add_subplot(gs[1, 2])]\n",
-    "xlims = [(0., 1.), (0., 1.), (0., 1.), (0., 1.)]\n",
+    "axes = [\n",
+    "    fig.add_subplot(gs[0, :]),\n",
+    "    fig.add_subplot(gs[1, 0]),\n",
+    "    fig.add_subplot(gs[1, 1]),\n",
+    "    fig.add_subplot(gs[1, 2]),\n",
+    "]\n",
+    "xlims = [(0.0, 1.0), (0.0, 1.0), (0.0, 1.0), (0.0, 1.0)]\n",
     "ylims = [(ymax, ymin), (s0, s1), (0, s1), (0, 600)]\n",
     "\n",
-    "letters = ['A', 'B', 'C', 'D']\n",
+    "letters = [\"A\", \"B\", \"C\", \"D\"]\n",
     "for idx, ax in enumerate(axes):\n",
     "    ax.set_xlim(xlims[idx])\n",
     "    ax.set_ylim(ylims[idx])\n",
     "    if idx > 0:\n",
     "        ax.set_xticks([0, 0.25, 0.5, 0.75, 1])\n",
-    "        ax.set_xticklabels(['0', '0.25', '0.50', '0.75', '1.00'])\n",
+    "        ax.set_xticklabels([\"0\", \"0.25\", \"0.50\", \"0.75\", \"1.00\"])\n",
     "    else:\n",
-    "        ax.tick_params(axis='both', length=0, labelbottom=False, labelleft=False) \n",
+    "        ax.tick_params(axis=\"both\", length=0, labelbottom=False, labelleft=False)\n",
     "    if spnspecs is not None:\n",
     "        spnspecs.remove_edge_ticks(ax)\n",
     "        spnspecs.heading(ax, letter=letters[idx])\n",
     "\n",
     "ax = axes[0]\n",
-    "ax.fill_between(dc, ymax, 0, color='0.9')\n",
+    "ax.fill_between(dc, ymax, 0, color=\"0.9\")\n",
     "# ax.plot(dc, z, lw=1, color='black')\n",
-    "ax.axhline(zmin, lw=0.5, ls='--', color='black')\n",
-    "ax.axhline(zmax, lw=.5, ls='--', color='black')\n",
-    "ax.axhline(0, lw=0.75, color='black')\n",
+    "ax.axhline(zmin, lw=0.5, ls=\"--\", color=\"black\")\n",
+    "ax.axhline(zmax, lw=0.5, ls=\"--\", color=\"black\")\n",
+    "ax.axhline(0, lw=0.75, color=\"black\")\n",
     "if spnspecs is not None:\n",
-    "    text = 'DDRN'\n",
-    "    #spnspecs.add_text(ax, 'DDRN', x=1.01, y=0.5, ha='center', va='center', \n",
+    "    text = \"DDRN\"\n",
+    "    # spnspecs.add_text(ax, 'DDRN', x=1.01, y=0.5, ha='center', va='center',\n",
     "    #                  bold=False, rotation=-90)\n",
-    "    spnspecs.add_annotation(ax, text=text, xy=(1.01, 0.77), xytext=(1.01, 0.5),\n",
-    "                            bold=False, rotation=-90, ha='center', va='center',\n",
-    "                            xycoords='axes fraction', textcoords='axes fraction',\n",
-    "                            arrowprops=arrowprops)\n",
-    "    spnspecs.add_annotation(ax, text=text, xy=(1.01, 0.23), xytext=(1.01, 0.5),\n",
-    "                            bold=False, rotation=-90, ha='center', va='center',\n",
-    "                            xycoords='axes fraction', textcoords='axes fraction',\n",
-    "                            arrowprops=arrowprops)\n",
-    "    text = 'Land surface elevation'\n",
-    "#     spnspecs.add_annotation(ax, text=text, xy=(0.6, 0.48), xytext=(0.68, 0.9),\n",
-    "#                             bold=False, ha='center', va='center',\n",
-    "#                             xycoords='axes fraction', textcoords='axes fraction',\n",
-    "#                             arrowprops=arrowprops)\n",
-    "    spnspecs.add_text(ax, text=text, x=0.5, y=0.5, bold=False, ha='center', va='bottom')\n",
-    "    text = 'Cell top'\n",
-    "    spnspecs.add_text(ax, text=text, x=0.01, y=0.5, bold=False, ha='left', va='bottom')\n",
-    "    text = 'Cell bottom'\n",
-    "    spnspecs.add_text(ax, text=text, x=0.01, y=0.01, bold=False, ha='left', va='bottom')\n",
-    "    text = r'Land surface elevation + $\\frac{DDRN}{2}$'\n",
-    "    spnspecs.add_text(ax, text=text, x=0.5, y=zmax, transform=False, bold=False, ha='center')\n",
-    "    text = r'Land surface elevation - $\\frac{DDRN}{2}$'\n",
-    "    spnspecs.add_text(ax, text=text, x=0.5, y=zmin, transform=False, bold=False, \n",
-    "                      ha='center', va='top')\n",
-    "    text = r'$ZDRN + DDRN$'\n",
-    "    spnspecs.add_text(ax, text=text, x=0.99, y=zmax, transform=False, bold=False, ha='right')\n",
-    "    text = r'$ZDRN$'\n",
-    "    spnspecs.add_text(ax, text=text, x=0.99, y=zmin, transform=False, bold=False, \n",
-    "                      ha='right', va='top')\n",
+    "    spnspecs.add_annotation(\n",
+    "        ax,\n",
+    "        text=text,\n",
+    "        xy=(1.01, 0.77),\n",
+    "        xytext=(1.01, 0.5),\n",
+    "        bold=False,\n",
+    "        rotation=-90,\n",
+    "        ha=\"center\",\n",
+    "        va=\"center\",\n",
+    "        xycoords=\"axes fraction\",\n",
+    "        textcoords=\"axes fraction\",\n",
+    "        arrowprops=arrowprops,\n",
+    "    )\n",
+    "    spnspecs.add_annotation(\n",
+    "        ax,\n",
+    "        text=text,\n",
+    "        xy=(1.01, 0.23),\n",
+    "        xytext=(1.01, 0.5),\n",
+    "        bold=False,\n",
+    "        rotation=-90,\n",
+    "        ha=\"center\",\n",
+    "        va=\"center\",\n",
+    "        xycoords=\"axes fraction\",\n",
+    "        textcoords=\"axes fraction\",\n",
+    "        arrowprops=arrowprops,\n",
+    "    )\n",
+    "    text = \"Land surface elevation\"\n",
+    "    #     spnspecs.add_annotation(ax, text=text, xy=(0.6, 0.48), xytext=(0.68, 0.9),\n",
+    "    #                             bold=False, ha='center', va='center',\n",
+    "    #                             xycoords='axes fraction', textcoords='axes fraction',\n",
+    "    #                             arrowprops=arrowprops)\n",
+    "    spnspecs.add_text(ax, text=text, x=0.5, y=0.5, bold=False, ha=\"center\", va=\"bottom\")\n",
+    "    text = \"Cell top\"\n",
+    "    spnspecs.add_text(ax, text=text, x=0.01, y=0.5, bold=False, ha=\"left\", va=\"bottom\")\n",
+    "    text = \"Cell bottom\"\n",
+    "    spnspecs.add_text(ax, text=text, x=0.01, y=0.01, bold=False, ha=\"left\", va=\"bottom\")\n",
+    "    text = r\"Land surface elevation + $\\frac{DDRN}{2}$\"\n",
+    "    spnspecs.add_text(\n",
+    "        ax, text=text, x=0.5, y=zmax, transform=False, bold=False, ha=\"center\"\n",
+    "    )\n",
+    "    text = r\"Land surface elevation - $\\frac{DDRN}{2}$\"\n",
+    "    spnspecs.add_text(\n",
+    "        ax, text=text, x=0.5, y=zmin, transform=False, bold=False, ha=\"center\", va=\"top\"\n",
+    "    )\n",
+    "    text = r\"$ZDRN + DDRN$\"\n",
+    "    spnspecs.add_text(\n",
+    "        ax, text=text, x=0.99, y=zmax, transform=False, bold=False, ha=\"right\"\n",
+    "    )\n",
+    "    text = r\"$ZDRN$\"\n",
+    "    spnspecs.add_text(\n",
+    "        ax, text=text, x=0.99, y=zmin, transform=False, bold=False, ha=\"right\", va=\"top\"\n",
+    "    )\n",
     "\n",
     "ax = axes[1]\n",
-    "ax.axhline(0, lw=0.5, ls='-.', color='black')\n",
-    "ax.axhline(1, lw=0.5, ls='-.', color='black')\n",
-    "ax.plot(t, s, lw=1.5, color='0.5', ls='--', label=r'$F_{DRN}$')\n",
-    "ax.set_xlabel('Fractional simulation time, unitless')\n",
-    "ax.set_ylabel(r'$\\frac{h - ZDRN}{DDRN}$, unitless')\n",
-    "text = r'$h = ZDRN + DDRN$'\n",
-    "spnspecs.add_text(ax, text=text, x=0.02, y=.98, transform=False, bold=False, ha='left', va='top')\n",
-    "text = r'$h = ZDRN$'\n",
-    "spnspecs.add_text(ax, text=text, x=0.99, y=0.01, transform=False, bold=False, \n",
-    "                  ha='right', va='bottom')\n",
+    "ax.axhline(0, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.axhline(1, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.plot(t, s, lw=1.5, color=\"0.5\", ls=\"--\", label=r\"$F_{DRN}$\")\n",
+    "ax.set_xlabel(\"Fractional simulation time, unitless\")\n",
+    "ax.set_ylabel(r\"$\\frac{h - ZDRN}{DDRN}$, unitless\")\n",
+    "text = r\"$h = ZDRN + DDRN$\"\n",
+    "spnspecs.add_text(\n",
+    "    ax, text=text, x=0.02, y=0.98, transform=False, bold=False, ha=\"left\", va=\"top\"\n",
+    ")\n",
+    "text = r\"$h = ZDRN$\"\n",
+    "spnspecs.add_text(\n",
+    "    ax, text=text, x=0.99, y=0.01, transform=False, bold=False, ha=\"right\", va=\"bottom\"\n",
+    ")\n",
     "\n",
     "\n",
     "ax = axes[2]\n",
-    "ax.plot(t, q0, lw=0.75, color='black')\n",
-    "ax.plot(t, ql, lw=1.5, color='blue')\n",
-    "ax.plot(t, qc, lw=0.75, color='red', label=r'$F_{DRN}^*$')\n",
-    "ax.set_xlabel('Fractional simulation time, unitless')\n",
-    "ax.set_ylabel(r'Drain discharge rate, , L$^3$/T')\n",
-    "text = r'Area = 1 L$^2$' + '\\n' + r'K$_v$ = 1 L/T'\n",
+    "ax.plot(t, q0, lw=0.75, color=\"black\")\n",
+    "ax.plot(t, ql, lw=1.5, color=\"blue\")\n",
+    "ax.plot(t, qc, lw=0.75, color=\"red\", label=r\"$F_{DRN}^*$\")\n",
+    "ax.set_xlabel(\"Fractional simulation time, unitless\")\n",
+    "ax.set_ylabel(r\"Drain discharge rate, , L$^3$/T\")\n",
+    "text = r\"Area = 1 L$^2$\" + \"\\n\" + r\"K$_v$ = 1 L/T\"\n",
     "if spnspecs is not None:\n",
-    "    spnspecs.add_text(ax, text=text, x=0.1, y=0.93, italic=False, bold=False, ha='left', va='top')\n",
+    "    spnspecs.add_text(\n",
+    "        ax, text=text, x=0.1, y=0.93, italic=False, bold=False, ha=\"left\", va=\"top\"\n",
+    "    )\n",
     "\n",
     "ax = axes[3]\n",
-    "ax.plot(t, q0.cumsum(), lw=0.75, color='black', label=r'$F_{DRN}^0$ (original)')\n",
-    "ax.plot(t, ql.cumsum(), lw=1.5, color='blue', label=r'$F_{DRN}$ (linear)')\n",
-    "ax.plot(t, qc.cumsum(), lw=0.75, color='red', label=r'$F_{DRN}^*$ (cubic)')\n",
-    "ax.set_xlabel('Fractional simulation time, unitless')\n",
-    "ax.set_ylabel(r'Cumulative drain discharge, L$^3$')\n",
+    "ax.plot(t, q0.cumsum(), lw=0.75, color=\"black\", label=r\"$F_{DRN}^0$ (original)\")\n",
+    "ax.plot(t, ql.cumsum(), lw=1.5, color=\"blue\", label=r\"$F_{DRN}$ (linear)\")\n",
+    "ax.plot(t, qc.cumsum(), lw=0.75, color=\"red\", label=r\"$F_{DRN}^*$ (cubic)\")\n",
+    "ax.set_xlabel(\"Fractional simulation time, unitless\")\n",
+    "ax.set_ylabel(r\"Cumulative drain discharge, L$^3$\")\n",
     "if spnspecs is not None:\n",
-    "    spnspecs.graph_legend(ax, loc='upper left', labelspacing=0.15)\n",
+    "    spnspecs.graph_legend(ax, loc=\"upper left\", labelspacing=0.15)\n",
     "\n",
-    "fpth = os.path.join(figpth, 'DRNDischargeDifferences.pdf')\n",
+    "fpth = os.path.join(figpth, \"DRNDischargeDifferences.pdf\")\n",
     "fig.savefig(fpth, dpi=dpi)"
    ]
   },
diff --git a/doc/SuppTechInfo/python/MAW-Discharge-Correction.ipynb b/doc/SuppTechInfo/python/MAW-Discharge-Correction.ipynb
index 4e02b435624..3520fcb8c27 100644
--- a/doc/SuppTechInfo/python/MAW-Discharge-Correction.ipynb
+++ b/doc/SuppTechInfo/python/MAW-Discharge-Correction.ipynb
@@ -7,9 +7,9 @@
    "outputs": [],
    "source": [
     "import os\n",
-    "import numpy as np\n",
+    "\n",
     "import matplotlib.pyplot as plt\n",
-    "import matplotlib.gridspec as gridspec"
+    "import numpy as np"
    ]
   },
   {
@@ -20,6 +20,7 @@
    "source": [
     "try:\n",
     "    import spnspecs\n",
+    "\n",
     "    spnspecs.set_graph_specifications()\n",
     "except:\n",
     "    spnspecs = None"
@@ -31,7 +32,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "figpth = '../Figures'\n",
+    "figpth = \"../Figures\"\n",
     "width = 6.8\n",
     "dpi = 300"
    ]
@@ -73,7 +74,7 @@
     "step = 0.01\n",
     "h0, h1 = -2, 2\n",
     "h = np.arange(h0, h1, step=step, dtype=dtype)\n",
-    "e = 0."
+    "e = 0.0"
    ]
   },
   {
@@ -91,7 +92,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hbar0 = drycell_smoothing(h, e, 0.)"
+    "hbar0 = drycell_smoothing(h, e, 0.0)"
    ]
   },
   {
@@ -100,7 +101,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hbar = drycell_smoothing(h, e, 2.)"
+    "hbar = drycell_smoothing(h, e, 2.0)"
    ]
   },
   {
@@ -109,7 +110,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(h,h)\n",
+    "plt.plot(h, h)\n",
     "plt.plot(h, hbar0)\n",
     "plt.plot(h, hbar)"
    ]
@@ -127,39 +128,58 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig, axes = plt.subplots(nrows=1, ncols=1, tight_layout=True, \n",
-    "                         figsize=(width/2, width/2))\n",
+    "fig, axes = plt.subplots(\n",
+    "    nrows=1, ncols=1, tight_layout=True, figsize=(width / 2, width / 2)\n",
+    ")\n",
     "if not isinstance(axes, list):\n",
     "    axes = [axes]\n",
     "\n",
-    "letters = ['A', 'B']\n",
+    "letters = [\"A\", \"B\"]\n",
     "for idx, ax in enumerate(axes):\n",
     "    ax.set_xlim(h0, h1)\n",
     "    ax.set_ylim(h0, h1)\n",
     "    ax.set_xticks([-1, 0, 1])\n",
-    "    ax.set_xticklabels([r'$e_{MAW,n} - \\epsilon$', r'$e_{MAW,n}$', r'$e_{MAW,n} + \\epsilon$'])\n",
+    "    ax.set_xticklabels(\n",
+    "        [r\"$e_{MAW,n} - \\epsilon$\", r\"$e_{MAW,n}$\", r\"$e_{MAW,n} + \\epsilon$\"]\n",
+    "    )\n",
     "    ax.set_yticks([0])\n",
-    "    ax.set_yticklabels([r'$e_{MAW,n}$'], rotation=90, va='center')\n",
+    "    ax.set_yticklabels([r\"$e_{MAW,n}$\"], rotation=90, va=\"center\")\n",
     "    if spnspecs is not None:\n",
     "        spnspecs.remove_edge_ticks(ax)\n",
-    "        #spnspecs.heading(ax, letter=letters[idx])\n",
-    "    \n",
+    "        # spnspecs.heading(ax, letter=letters[idx])\n",
+    "\n",
     "ax = axes[0]\n",
-    "ax.set_aspect('equal', 'box')\n",
-    "ax.axhline(0, lw=0.5, ls='-.', color='black')\n",
-    "ax.axvline(0, lw=0.5, ls='-.', color='black')\n",
-    "ax.plot([-3], [-3], lw=0, mfc='black', mec='none', ms=4.5, marker=r'$\\epsilon$', label='smothing distance')\n",
-    "ax.plot(h, hbar, lw=1.25, color='#00BFFF', label=r'$\\overline{h}_n^*$ continuous')\n",
-    "ax.plot(h, hbar0, lw=1., color='black', ls='-', label=r'$\\overline{h}_n$ discontinuous')\n",
+    "ax.set_aspect(\"equal\", \"box\")\n",
+    "ax.axhline(0, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.axvline(0, lw=0.5, ls=\"-.\", color=\"black\")\n",
+    "ax.plot(\n",
+    "    [-3],\n",
+    "    [-3],\n",
+    "    lw=0,\n",
+    "    mfc=\"black\",\n",
+    "    mec=\"none\",\n",
+    "    ms=4.5,\n",
+    "    marker=r\"$\\epsilon$\",\n",
+    "    label=\"smothing distance\",\n",
+    ")\n",
+    "ax.plot(h, hbar, lw=1.25, color=\"#00BFFF\", label=r\"$\\overline{h}_n^*$ continuous\")\n",
+    "ax.plot(\n",
+    "    h, hbar0, lw=1.0, color=\"black\", ls=\"-\", label=r\"$\\overline{h}_n$ discontinuous\"\n",
+    ")\n",
     "if spnspecs is not None:\n",
     "    handles, labels = ax.get_legend_handles_labels()\n",
-    "    spnspecs.graph_legend(ax, handles=handles[::-1], labels=labels[::-1],\n",
-    "                          loc='lower right', bbox_to_anchor=(1.00,0.05),\n",
-    "                          handlelength=1.5)\n",
-    "ax.set_xlabel('Downstream head, in length units')\n",
-    "ax.set_ylabel('Reference elevation, in length units')\n",
+    "    spnspecs.graph_legend(\n",
+    "        ax,\n",
+    "        handles=handles[::-1],\n",
+    "        labels=labels[::-1],\n",
+    "        loc=\"lower right\",\n",
+    "        bbox_to_anchor=(1.00, 0.05),\n",
+    "        handlelength=1.5,\n",
+    "    )\n",
+    "ax.set_xlabel(\"Downstream head, in length units\")\n",
+    "ax.set_ylabel(\"Reference elevation, in length units\")\n",
     "\n",
-    "fpth = os.path.join(figpth, 'MAWDischargeCorrection.pdf')\n",
+    "fpth = os.path.join(figpth, \"MAWDischargeCorrection.pdf\")\n",
     "fig.savefig(fpth, dpi=dpi)"
    ]
   },
diff --git a/doc/SuppTechInfo/python/STO-SpecificStorage-Derivatives.ipynb b/doc/SuppTechInfo/python/STO-SpecificStorage-Derivatives.ipynb
index 8da9aa19418..bbbbf7e82b1 100644
--- a/doc/SuppTechInfo/python/STO-SpecificStorage-Derivatives.ipynb
+++ b/doc/SuppTechInfo/python/STO-SpecificStorage-Derivatives.ipynb
@@ -6,8 +6,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import numpy as np\n",
-    "import matplotlib.pyplot as plt"
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np"
    ]
   },
   {
@@ -17,9 +17,9 @@
    "outputs": [],
    "source": [
     "nz = 10000\n",
-    "cell_size = 10.\n",
+    "cell_size = 10.0\n",
     "sat = 0.25\n",
-    "bot = 10.\n",
+    "bot = 10.0\n",
     "h = bot + sat * cell_size\n",
     "top = bot + cell_size\n",
     "deltaz = top - bot"
@@ -31,7 +31,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "z = np.linspace(bot, h, nz+1)"
+    "z = np.linspace(bot, h, nz + 1)"
    ]
   },
   {
@@ -40,7 +40,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ss = 1."
+    "ss = 1.0"
    ]
   },
   {
@@ -49,9 +49,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ss_fd = 0.\n",
+    "ss_fd = 0.0\n",
     "for idx in range(nz):\n",
-    "    z0, z1 = z[idx], z[idx + 1] \n",
+    "    z0, z1 = z[idx], z[idx + 1]\n",
     "    dz = z1 - z0\n",
     "    rho = ss * dz\n",
     "    psi = h - z0 - 0.5 * dz\n",
@@ -65,7 +65,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "psi_avg = (h - bot - deltaz * sat / 2)\n",
+    "psi_avg = h - bot - deltaz * sat / 2\n",
     "psi_avg"
    ]
   },
@@ -98,26 +98,26 @@
     "        x = np.array([x], dtype=float)\n",
     "    elif isinstance(x, (list, tuple)):\n",
     "        x = np.array(x, dtype=float)\n",
-    "        \n",
+    "\n",
     "    y = np.zeros(x.shape, dtype=z.dtype)\n",
     "    dz = tp - bt\n",
     "    br = (x - bt) / dz\n",
-    "    br[x < bt] = 0.\n",
-    "    br[x > tp] = 1.\n",
-    "    \n",
+    "    br[x < bt] = 0.0\n",
+    "    br[x > tp] = 1.0\n",
+    "\n",
     "    if omega == 0:\n",
     "        y[:] = br\n",
     "    else:\n",
-    "        av = 1. / (1. - omega)\n",
+    "        av = 1.0 / (1.0 - omega)\n",
     "\n",
     "        idx = br < omega\n",
-    "        y[idx] = 0.5 * av * br[idx]**2. / omega \n",
+    "        y[idx] = 0.5 * av * br[idx] ** 2.0 / omega\n",
     "\n",
-    "        idx = (br >= omega) & (br < 1. - omega)\n",
-    "        y[idx] = av * br[idx] + 0.5 * (1. - av)\n",
+    "        idx = (br >= omega) & (br < 1.0 - omega)\n",
+    "        y[idx] = av * br[idx] + 0.5 * (1.0 - av)\n",
     "\n",
-    "        idx = (br >= 1. - omega)\n",
-    "        y[idx] = 1. - ((0.5 * av * (1. - br[idx])**2) / omega)\n",
+    "        idx = br >= 1.0 - omega\n",
+    "        y[idx] = 1.0 - ((0.5 * av * (1.0 - br[idx]) ** 2) / omega)\n",
     "\n",
     "    return y"
    ]
@@ -132,12 +132,11 @@
     "    dx = 1e-6\n",
     "    derv = np.zeros(x.shape, dtype=x.dtype)\n",
     "    for idx, xx in enumerate(x):\n",
-    "        xx0, xx1 = xx - dx, xx + dx \n",
+    "        xx0, xx1 = xx - dx, xx + dx\n",
     "        y0 = quad_sat(xx0, tp, bt, omega=omega)\n",
     "        y1 = quad_sat(xx1, tp, bt, omega=omega)\n",
     "        derv[idx] = (y1**power - y0**power) / (2 * dx)\n",
-    "    return derv\n",
-    "        "
+    "    return derv"
    ]
   },
   {
@@ -146,7 +145,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "harr = np.linspace(bot - 1., top + 1., nz)\n",
+    "harr = np.linspace(bot - 1.0, top + 1.0, nz)\n",
     "omega = 1e-6"
    ]
   },
@@ -156,7 +155,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sat_lin = quad_sat(harr, top, bot, omega=0.)\n",
+    "sat_lin = quad_sat(harr, top, bot, omega=0.0)\n",
     "sat_lin.mean()"
    ]
   },
diff --git a/doc/SuppTechInfo/python/STO-SpecificStorage.ipynb b/doc/SuppTechInfo/python/STO-SpecificStorage.ipynb
index c01970c492b..90259af5d70 100644
--- a/doc/SuppTechInfo/python/STO-SpecificStorage.ipynb
+++ b/doc/SuppTechInfo/python/STO-SpecificStorage.ipynb
@@ -7,8 +7,9 @@
    "outputs": [],
    "source": [
     "import os\n",
-    "import numpy as np\n",
-    "import matplotlib.pyplot as plt"
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np"
    ]
   },
   {
@@ -18,6 +19,7 @@
    "outputs": [],
    "source": [
     "import spnspecs\n",
+    "\n",
     "spnspecs.set_graph_specifications()"
    ]
   },
@@ -27,7 +29,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "figpth = '../Figures'\n",
+    "figpth = \"../Figures\"\n",
     "width = 6.8\n",
     "dpi = 300"
    ]
@@ -38,7 +40,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#spnspecs.set_graph_specifications()"
+    "# spnspecs.set_graph_specifications()"
    ]
   },
   {
@@ -47,7 +49,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def calc_err(tp, datum=1000., ss=1e-5, bt=0., sf0=0.01, sf1=0.99, dt=1., area=1.):\n",
+    "def calc_err(tp, datum=1000.0, ss=1e-5, bt=0.0, sf0=0.01, sf1=0.99, dt=1.0, area=1.0):\n",
     "    dz = tp - bt\n",
     "    h0 = sf0 * dz + bt\n",
     "    h1 = sf1 * dz + bt\n",
@@ -65,7 +67,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def err_function(tp, datum=1000., ss=1e-5, bt=0., sf0=0.01, sf1=0.99, dt=1., area=1.):\n",
+    "def err_function(\n",
+    "    tp, datum=1000.0, ss=1e-5, bt=0.0, sf0=0.01, sf1=0.99, dt=1.0, area=1.0\n",
+    "):\n",
     "    dz = tp - bt\n",
     "    rho1 = ss * area * dz / dt\n",
     "    bt1 = bt + datum\n",
@@ -78,7 +82,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def err_function_simp(tp, datum=1000., ss=1e-5, bt=0., sf0=0.01, sf1=0.99, dt=1., area=1.):\n",
+    "def err_function_simp(\n",
+    "    tp, datum=1000.0, ss=1e-5, bt=0.0, sf0=0.01, sf1=0.99, dt=1.0, area=1.0\n",
+    "):\n",
     "    dz = tp - bt\n",
     "    rho1 = ss * area * dz / dt\n",
     "    return -rho1 * datum * (sf0 - sf1)"
@@ -90,7 +96,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "top = np.array([0.1, 0.5, 1., 5., 10., 50., 100., 500., 1000.])\n",
+    "top = np.array([0.1, 0.5, 1.0, 5.0, 10.0, 50.0, 100.0, 500.0, 1000.0])\n",
     "err = calc_err(top, ss=1e-5)"
    ]
   },
@@ -110,24 +116,23 @@
    "outputs": [],
    "source": [
     "fig, axes = plt.subplots(\n",
-    "    nrows=1, \n",
-    "    ncols=1, \n",
-    "    constrained_layout=True, \n",
+    "    nrows=1,\n",
+    "    ncols=1,\n",
+    "    constrained_layout=True,\n",
     "    figsize=(width, 0.3 * width),\n",
     ")\n",
     "\n",
     "ax = axes\n",
-    "for (sf0, color) in zip(\n",
-    "    (0.0, 0.1, 0.2, 0.3, 0.4, 0.5), \n",
-    "    (\"red\", \"orange\", \"green\", \"cyan\", \"blue\", \"black\")\n",
+    "for sf0, color in zip(\n",
+    "    (0.0, 0.1, 0.2, 0.3, 0.4, 0.5), (\"red\", \"orange\", \"green\", \"cyan\", \"blue\", \"black\")\n",
     "):\n",
-    "    sf1 = 1. - sf0\n",
+    "    sf1 = 1.0 - sf0\n",
     "    label = r\"S$_F^t$ = {:.2f} S$_F^{{told}}$ = {:.2f}\".format(sf1, sf0)\n",
     "    ax.plot(\n",
-    "        top, \n",
-    "        calc_err(top, ss=1, sf0=sf0, sf1=sf1), \n",
-    "        lw=1.25, \n",
-    "        color=color, \n",
+    "        top,\n",
+    "        calc_err(top, ss=1, sf0=sf0, sf1=sf1),\n",
+    "        lw=1.25,\n",
+    "        color=color,\n",
     "        label=label,\n",
     "        clip_on=False,\n",
     "    )\n",
@@ -136,7 +141,7 @@
     "ax.set_xlabel(\"Cell thickness, in length units\")\n",
     "ax.set_ylabel(r\"Q$_{SS}$\" + \" Differences,\\nin cubic length per time units\")\n",
     "spnspecs.graph_legend(\n",
-    "    ax, \n",
+    "    ax,\n",
     "    ncol=2,\n",
     "    labelspacing=0.2,\n",
     "    columnspacing=0.5,\n",
@@ -153,21 +158,21 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig, axes = plt.subplots(nrows=1, ncols=1, tight_layout=True, \n",
-    "                         figsize=(width, width/2))\n",
+    "fig, axes = plt.subplots(\n",
+    "    nrows=1, ncols=1, tight_layout=True, figsize=(width, width / 2)\n",
+    ")\n",
     "\n",
     "ax = axes\n",
-    "for (sf0, color) in zip(\n",
-    "    (0.0, 0.1, 0.2, 0.3, 0.4, 0.5), \n",
-    "    (\"red\", \"orange\", \"green\", \"cyan\", \"blue\", \"black\")\n",
+    "for sf0, color in zip(\n",
+    "    (0.0, 0.1, 0.2, 0.3, 0.4, 0.5), (\"red\", \"orange\", \"green\", \"cyan\", \"blue\", \"black\")\n",
     "):\n",
-    "    sf1 = 1. - sf0\n",
+    "    sf1 = 1.0 - sf0\n",
     "    label = r\"S$_F^t$ = {:.2f} S$_F^{{told}}$ = {:.2f}\".format(sf1, sf0)\n",
     "    ax.plot(\n",
-    "        top, \n",
-    "        err_function(top, ss=1, sf0=sf0, sf1=sf1), \n",
-    "        lw=1.25, \n",
-    "        color=color, \n",
+    "        top,\n",
+    "        err_function(top, ss=1, sf0=sf0, sf1=sf1),\n",
+    "        lw=1.25,\n",
+    "        color=color,\n",
     "        label=label,\n",
     "        clip_on=False,\n",
     "    )\n",
@@ -186,24 +191,23 @@
    "outputs": [],
    "source": [
     "fig, axes = plt.subplots(\n",
-    "    nrows=1, \n",
-    "    ncols=1, \n",
-    "    constrained_layout=True, \n",
+    "    nrows=1,\n",
+    "    ncols=1,\n",
+    "    constrained_layout=True,\n",
     "    figsize=(width, 0.3 * width),\n",
     ")\n",
     "\n",
     "ax = axes\n",
-    "for (sf0, color) in zip(\n",
-    "    (0.0, 0.1, 0.2, 0.3, 0.4, 0.5), \n",
-    "    (\"red\", \"orange\", \"green\", \"cyan\", \"blue\", \"black\")\n",
+    "for sf0, color in zip(\n",
+    "    (0.0, 0.1, 0.2, 0.3, 0.4, 0.5), (\"red\", \"orange\", \"green\", \"cyan\", \"blue\", \"black\")\n",
     "):\n",
-    "    sf1 = 1. - sf0\n",
+    "    sf1 = 1.0 - sf0\n",
     "    label = r\"S$_F^t$ = {:.2f} S$_F^{{told}}$ = {:.2f}\".format(sf1, sf0)\n",
     "    ax.plot(\n",
-    "        top, \n",
-    "        err_function_simp(top, ss=1e-5, sf0=sf0, sf1=sf1), \n",
-    "        lw=1.25, \n",
-    "        color=color, \n",
+    "        top,\n",
+    "        err_function_simp(top, ss=1e-5, sf0=sf0, sf1=sf1),\n",
+    "        lw=1.25,\n",
+    "        color=color,\n",
     "        label=label,\n",
     "        clip_on=False,\n",
     "    )\n",
@@ -212,14 +216,13 @@
     "ax.set_xlabel(\"Cell thickness, in meters\")\n",
     "ax.set_ylabel(r\"$\\Delta$ Q$_{SS}$, in meters$^3$/ day\")\n",
     "spnspecs.graph_legend(\n",
-    "    ax, \n",
+    "    ax,\n",
     "    ncol=2,\n",
     "    labelspacing=0.2,\n",
     "    columnspacing=0.5,\n",
     ")\n",
-    "spnspecs.remove_edge_ticks(ax);\n",
-    "\n",
-    "fpth = os.path.join(figpth, 'STOSsError.pdf')\n",
+    "spnspecs.remove_edge_ticks(ax)\n",
+    "fpth = os.path.join(figpth, \"STOSsError.pdf\")\n",
     "fig.savefig(fpth, dpi=dpi);"
    ]
   },
diff --git a/doc/SuppTechInfo/python/VSC-nonlinear.ipynb b/doc/SuppTechInfo/python/VSC-nonlinear.ipynb
index cd69f368048..21e897dc32b 100644
--- a/doc/SuppTechInfo/python/VSC-nonlinear.ipynb
+++ b/doc/SuppTechInfo/python/VSC-nonlinear.ipynb
@@ -6,11 +6,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import os\n",
-    "import numpy as np\n",
+    "import math\n",
+    "\n",
     "import matplotlib.pyplot as plt\n",
-    "import matplotlib.gridspec as gridspec\n",
-    "import math"
+    "import numpy as np"
    ]
   },
   {
@@ -21,6 +20,7 @@
    "source": [
     "try:\n",
     "    import spnspecs\n",
+    "\n",
     "    spnspecs.set_graph_specifications()\n",
     "except:\n",
     "    spnspecs = None"
@@ -32,7 +32,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "figpth = '../Figures'\n",
+    "figpth = \"../Figures\"\n",
     "width = 6.8\n",
     "dpi = 300"
    ]
@@ -73,7 +73,7 @@
     "    for idx, val in enumerate(T):\n",
     "        mu_T = a1 * a2 ** (a3 / (T[idx] + a4))\n",
     "        visc[idx] = mu_T\n",
-    "        \n",
+    "\n",
     "    return visc"
    ]
   },
@@ -104,13 +104,14 @@
     "# reference viscosity\n",
     "mu_0 = nonlin_vsc(np.array([Tviscref]))\n",
     "\n",
+    "\n",
     "def lin_vsc(T, dviscdT, Tviscref):\n",
     "    # use the linear formula\n",
     "    visc_lin = np.zeros(T.shape, dtype=np.float)\n",
     "    for idx, val in enumerate(T):\n",
     "        mu_T = dviscdT * (val - Tviscref) + mu_0\n",
     "        visc_lin[idx] = mu_T\n",
-    "    \n",
+    "\n",
     "    return visc_lin"
    ]
   },
@@ -130,8 +131,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(T,visc)\n",
-    "plt.plot(T,visc_lin)"
+    "plt.plot(T, visc)\n",
+    "plt.plot(T, visc_lin)"
    ]
   },
   {
@@ -150,7 +151,11 @@
     "x_increment = 5\n",
     "x_ticks = np.arange(T0, T1 + x_increment, x_increment)\n",
     "y_increment = 2e-4\n",
-    "y_ticks = np.arange(math.floor(visc.min() * 10000) / 10000, math.ceil(visc.max() * 10000) / 10000 + y_increment, y_increment)\n",
+    "y_ticks = np.arange(\n",
+    "    math.floor(visc.min() * 10000) / 10000,\n",
+    "    math.ceil(visc.max() * 10000) / 10000 + y_increment,\n",
+    "    y_increment,\n",
+    ")\n",
     "print(x_ticks)\n",
     "print(y_ticks)"
    ]
@@ -161,13 +166,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#fig, axes = plt.subplots(nrows=1, ncols=1, tight_layout=True, \n",
+    "# fig, axes = plt.subplots(nrows=1, ncols=1, tight_layout=True,\n",
     "#                       figsize=(width/1.5, width/2))\n",
     "\n",
-    "#if not isinstance(axes, list):\n",
+    "# if not isinstance(axes, list):\n",
     "#    axes = [axes]\n",
     "\n",
-    "#for idx, ax in enumerate(axes):\n",
+    "# for idx, ax in enumerate(axes):\n",
     "#    ax.set_xlim(T0, T1)\n",
     "#    ax.set_ylim(y_ticks[0], y_ticks[-2])\n",
     "#    ax.set_xticks(x_ticks)\n",
@@ -177,30 +182,30 @@
     "#    if spnspecs is not None:\n",
     "#        spnspecs.remove_edge_ticks(ax)\n",
     "#        spnspecs.heading(ax, letter=letters[idx])\n",
-    "    \n",
     "\n",
-    "#ax = axes[0]\n",
+    "\n",
+    "# ax = axes[0]\n",
     "##ax.set_aspect('square', adjustable='box')\n",
     "##ax.axhline(0, lw=0.5, ls='-.', color='black')\n",
     "##ax.axvline(0, lw=0.5, ls='-.', color='black')\n",
     "\n",
-    "#ax.plot(T, visc, lw=1.25, color='#00BFFF', label=r'$\\mu')\n",
-    "#ax.plot(T, visc_lin, lw=1., color='black', ls='-', label=r'$\\frac {\\partial \\mu} {\\partial T}')\n",
-    "#ax.plot([Tviscref], [mu_0], marker=\"o\", markersize=4, markeredgecolor=\"red\", markerfacecolor=\"blue\")\n",
-    "#text = r'$\\left ( T_0, \\mu_0 \\right )$'\n",
+    "# ax.plot(T, visc, lw=1.25, color='#00BFFF', label=r'$\\mu')\n",
+    "# ax.plot(T, visc_lin, lw=1., color='black', ls='-', label=r'$\\frac {\\partial \\mu} {\\partial T}')\n",
+    "# ax.plot([Tviscref], [mu_0], marker=\"o\", markersize=4, markeredgecolor=\"red\", markerfacecolor=\"blue\")\n",
+    "# text = r'$\\left ( T_0, \\mu_0 \\right )$'\n",
     "##ax.text(21, 0.00102, text, style='italic') #fontsize=6)\n",
-    "#spnspecs.add_text(ax, text=text, x=0.5, y=0.5, transform=False, bold=False, ha='left', va='top')\n",
+    "# spnspecs.add_text(ax, text=text, x=0.5, y=0.5, transform=False, bold=False, ha='left', va='top')\n",
     "\n",
-    "#if spnspecs is not None:\n",
+    "# if spnspecs is not None:\n",
     "#    handles, labels = ax.get_legend_handles_labels()\n",
     "#    spnspecs.graph_legend(ax, handles=handles[::-1], labels=labels[::-1],\n",
     "#                          loc='upper right', bbox_to_anchor=(1.00,0.95),\n",
     "#                          handlelength=1.5)\n",
-    "#ax.set_xlabel(r'Temperature, in $^\\circ$C')\n",
-    "#ax.set_ylabel(r'$\\mu$, in $\\frac{kg}{m \\cdot s}$')\n",
+    "# ax.set_xlabel(r'Temperature, in $^\\circ$C')\n",
+    "# ax.set_ylabel(r'$\\mu$, in $\\frac{kg}{m \\cdot s}$')\n",
     "\n",
-    "#fpth = os.path.join(figpth, 'VSCnonlinear.pdf')\n",
-    "#fig.savefig(fpth, dpi=dpi)"
+    "# fpth = os.path.join(figpth, 'VSCnonlinear.pdf')\n",
+    "# fig.savefig(fpth, dpi=dpi)"
    ]
   },
   {
diff --git a/ruff.toml b/ruff.toml
index 0cd8006214b..a8802328d7e 100644
--- a/ruff.toml
+++ b/ruff.toml
@@ -1,14 +1,5 @@
 line-length = 88
-target-version = "py37"
-include = [
-    ".doc/**/*.py",
-    "autotest/pyproject.toml",
-    "autotest/**/*.py",
-    "distribution/**/*.py",
-    "doc/**/*.py",
-    "pymake/**/*.py",
-    "utils/**/*.py",
-]
+target-version = "py39"
 
 [lint]
 select = [
@@ -24,6 +15,9 @@ ignore = [
     "F841", # local variable assigned but never used
 ]
 
+[lint.per-file-ignores]
+"doc/**.ipynb" = ["E501"]
+
 [format]
 quote-style = "double"
 indent-style = "space"