diff --git a/.build_rtd_docs/conf.py b/.build_rtd_docs/conf.py
index 1d4829f789e..1fbbefc769d 100644
--- a/.build_rtd_docs/conf.py
+++ b/.build_rtd_docs/conf.py
@@ -24,12 +24,12 @@
 on_rtd = os.environ.get("READTHEDOCS") == "True"
 
 # -- print current directory
-print("Current Directory...'{}'".format(os.path.abspath(os.getcwd())))
+print(f"Current Directory...'{os.path.abspath(os.getcwd())}'")
 
 # -- clean up doxygen files -------------------------------------------------
 dox_pths = ("_mf6io",)
 for dox_pth in dox_pths:
-    print("cleaning....{}".format(dox_pth))
+    print(f"cleaning....{dox_pth}")
     for root, dirs, files in os.walk(dox_pth):
         for name in files:
             fpth = os.path.join(root, name)
@@ -99,7 +99,8 @@
 if stdout:
     print(stdout.decode("utf-8"))
 if stderr:
-    print("Errors:\n{}".format(stderr.decode("utf-8")))
+    print("Errors:")
+    print(stderr.decode("utf-8"))
 
 # -- copy deprecations markdown ---------------------------------------------
 print("Copy the deprecations table")
@@ -120,7 +121,8 @@
 if stdout:
     print(stdout.decode("utf-8"))
 if stderr:
-    print("Errors:\n{}".format(stderr.decode("utf-8")))
+    print("Errors:")
+    print(stderr.decode("utf-8"))
 
 # -- update the doxygen version number ---------------------------------------
 print("Update the Doxyfile with the latest version number")
@@ -131,7 +133,7 @@
 with open("Doxyfile", "w") as fp:
     for line in lines:
         if tag in line:
-            line = '{}         = "version {}"\n'.format(tag, __version__)
+            line = f'{tag}         = "version {__version__}"\n'
         fp.write(line)
 
 # -- Project information -----------------------------------------------------
diff --git a/.doc/conf.py b/.doc/conf.py
index 26a18ee9c3b..79e6f19bb41 100644
--- a/.doc/conf.py
+++ b/.doc/conf.py
@@ -29,9 +29,9 @@
         src = os.path.join(src_pth, on_dir)
         dst = os.path.join(".", on_dir)
         if os.path.exists(dst):
-            print("deleting...{}".format(dst))
+            print(f"deleting...{dst}")
             shutil.rmtree(dst)
-        print("copying {} -> {}".format(src, dst))
+        print(f"copying {src} -> {dst}")
         shutil.copytree(src, dst)
 
     # copy files
@@ -40,9 +40,9 @@
         src = os.path.join(src_pth, file_name)
         dst = os.path.join(".", file_name)
         if os.path.exists(dst):
-            print("deleting...{}".format(dst))
+            print(f"deleting...{dst}")
             os.remove(dst)
-        print("copying {} -> {}".format(src, dst))
+        print(f"copying {src} -> {dst}")
         shutil.copy(src, dst)
 
 
diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs
index f2db7326990..c38610503fb 100644
--- a/.git-blame-ignore-revs
+++ b/.git-blame-ignore-revs
@@ -36,3 +36,6 @@ d56a12e0db8ff9d052b7ef9d7157506528e1a70e
 
 # reformat Python code with line length = 88 (#2056)
 f6c3a55ea72065ba1ce8540d1fd3c051696df1f1
+
+# reformat multi-line statements (#2070)
+cc7f930b961345536529942fc57db0c515a1f7fc
diff --git a/.github/common/update_compat_tables.py b/.github/common/update_compat_tables.py
index 59255741806..e5aab0f1bad 100644
--- a/.github/common/update_compat_tables.py
+++ b/.github/common/update_compat_tables.py
@@ -27,7 +27,7 @@
     ct = (
         "<!-- "
         + name
-        + " compat starts -->{}<!-- ".format("\n{}\n".format(table))
+        + " compat starts -->{}<!-- ".format(f"\n{table}\n")
         + name
         + " compat ends -->"
     )
diff --git a/.hpc/sstat_poll.py b/.hpc/sstat_poll.py
index def27f6ce7c..157c8c44048 100644
--- a/.hpc/sstat_poll.py
+++ b/.hpc/sstat_poll.py
@@ -1,6 +1,5 @@
 import argparse
 import pathlib as pl
-import sys
 import time
 from subprocess import PIPE, STDOUT, Popen
 
@@ -107,9 +106,6 @@ def _run_command(
     )
     slurm_args = parser.parse_args()
 
-    if sys.version_info < (3, 8):
-        sys.exit("Python version must be 3.8 or higher.")
-
     print(f"SLURM command: {slurm_args.command}")
     print(f"JobID: {slurm_args.jobid}")
 
diff --git a/autotest/common_regression.py b/autotest/common_regression.py
index 8114cc37e97..62bad292950 100644
--- a/autotest/common_regression.py
+++ b/autotest/common_regression.py
@@ -1,7 +1,8 @@
 import os
 import shutil
+from collections.abc import Iterator
 from pathlib import Path
-from typing import Iterator, List, Optional, Tuple, Union
+from typing import Optional, Union
 from warnings import warn
 
 COMPARE_PROGRAMS = (
@@ -221,8 +222,7 @@ def get_matching_files(
         extensions = [extensions]
 
     for ext in extensions:
-        for file in workspace.glob(f"*.{ext}"):
-            yield file
+        yield from workspace.glob(f"*.{ext}")
 
 
 def get_mf6_comparison(src):
@@ -475,7 +475,7 @@ def get_mf6_ftypes(namefile, ftypekeys):
 
 def get_regression_files(
     workspace: os.PathLike, extensions
-) -> Tuple[List[str], List[str]]:
+) -> tuple[list[str], list[str]]:
     if isinstance(extensions, str):
         extensions = [extensions]
     files = os.listdir(workspace)
diff --git a/autotest/conftest.py b/autotest/conftest.py
index 1d56dc869ad..c09eddcf199 100644
--- a/autotest/conftest.py
+++ b/autotest/conftest.py
@@ -1,6 +1,5 @@
 import sys
 from pathlib import Path
-from typing import Dict
 from warnings import warn
 
 import pytest
@@ -54,7 +53,7 @@ def bin_path() -> Path:
 
 
 @pytest.fixture(scope="session")
-def targets() -> Dict[str, Path]:
+def targets() -> dict[str, Path]:
     """
     Target executables for tests. These include local development builds as
     well as binaries 1) downloaded from GitHub and 2) rebuilt from the last
@@ -81,7 +80,7 @@ def targets() -> Dict[str, Path]:
     return d
 
 
-def try_get_target(targets: Dict[str, Path], name: str) -> Path:
+def try_get_target(targets: dict[str, Path], name: str) -> Path:
     """Try to retrieve the path to a binary. If the binary is a development
     target and can't be found, an error is raised. Otherwise (if the binary
     is downloaded or rebuilt) the test is skipped. This is to allow testing
diff --git a/autotest/framework.py b/autotest/framework.py
index c3e2ce32285..773394c7a31 100644
--- a/autotest/framework.py
+++ b/autotest/framework.py
@@ -1,10 +1,11 @@
 import os
 import shutil
+from collections.abc import Iterable
 from itertools import repeat
 from pathlib import Path
 from subprocess import PIPE, STDOUT, Popen
 from traceback import format_exc
-from typing import Callable, Dict, Iterable, List, Optional, Tuple, Union
+from typing import Callable, Optional, Union
 from warnings import warn
 
 import flopy
@@ -44,7 +45,7 @@
 )
 
 
-def api_return(success, model_ws) -> Tuple[bool, List[str]]:
+def api_return(success, model_ws) -> tuple[bool, list[str]]:
     """
     parse libmf6 stdout shared object file
     """
@@ -74,7 +75,7 @@ def get_workspace(sim_or_model) -> Path:
         raise ValueError(f"Unsupported model type: {type(sim_or_model)}")
 
 
-def run_parallel(workspace, target, ncpus) -> Tuple[bool, List[str]]:
+def run_parallel(workspace, target, ncpus) -> tuple[bool, list[str]]:
     if not is_in_ci() and get_ostag() in ["mac"]:
         oversubscribed = ["--hostfile", "localhost"]
         with open(f"{workspace}/localhost", "w") as f:
@@ -219,7 +220,7 @@ def __init__(
         self,
         name: str,
         workspace: Union[str, os.PathLike],
-        targets: Dict[str, Path],
+        targets: dict[str, Path],
         api_func: Optional[Callable] = None,
         build: Optional[Callable] = None,
         check: Optional[Callable] = None,
@@ -358,10 +359,8 @@ def _compare_heads(
                 verbose=self.verbose,
             )
             print(
-                (
-                    f"{EXTTEXT[extension]} comparison {i + 1}"
-                    + f"{self.name} ({os.path.basename(fpth0)})"
-                )
+                f"{EXTTEXT[extension]} comparison {i + 1}"
+                + f"{self.name} ({os.path.basename(fpth0)})"
             )
             if not success:
                 return False
@@ -540,7 +539,7 @@ def _run_sim_or_model(
         target: Union[str, os.PathLike],
         xfail: bool = False,
         ncpus: int = 1,
-    ) -> Tuple[bool, List[str]]:
+    ) -> tuple[bool, list[str]]:
         """
         Run a simulation or model with FloPy.
 
diff --git a/autotest/test_gwf_ats_lak01.py b/autotest/test_gwf_ats_lak01.py
index a7fecc5fd0c..3250acec556 100644
--- a/autotest/test_gwf_ats_lak01.py
+++ b/autotest/test_gwf_ats_lak01.py
@@ -400,9 +400,7 @@ def check_output(idx, test):
                     all_passed = False
                     msg = (
                         "recharge must be zero if overlying lake is "
-                        "active. node {} qlak {} qrch {} time {}".format(
-                            n0, qlakleak[n0], q, t
-                        )
+                        f"active. node {n0} qlak {qlakleak[n0]} qrch {q} time {t}"
                     )
                     print(msg)
     assert all_passed, "found recharge applied to cell beneath active lake"
diff --git a/autotest/test_gwf_ifmod_buy.py b/autotest/test_gwf_ifmod_buy.py
index 2e816e3c107..8a0b93707fb 100644
--- a/autotest/test_gwf_ifmod_buy.py
+++ b/autotest/test_gwf_ifmod_buy.py
@@ -555,52 +555,38 @@ def check_output(idx, test):
 
     # compare heads
     maxdiff = np.amax(abs(heads - heads_2models))
-    assert maxdiff < 10 * hclose_check, "Max. head diff. {} should \
-                     be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. head diff. {maxdiff} should \
+                     be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_x left
     maxdiff = np.amax(abs(qxb[:, :, 0:5] - qxb_left))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (x) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (x) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_y left
     maxdiff = np.amax(abs(qyb[:, :, 0:5] - qyb_left))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (y) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (y) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_z left
     maxdiff = np.amax(abs(qzb[:, :, 0:5] - qzb_left))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (z) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (z) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_x right
     maxdiff = np.amax(abs(qxb[:, :, 5:] - qxb_right))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (x) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (x) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_y right
     maxdiff = np.amax(abs(qyb[:, :, 5:] - qyb_right))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (y) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (y) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_z right
     maxdiff = np.amax(abs(qzb[:, :, 5:] - qzb_right))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (z) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (z) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # check budget error from .lst file
     for mname in [mname_ref, mname_left, mname_right]:
@@ -608,10 +594,9 @@ def check_output(idx, test):
         for line in open(fpth):
             if line.lstrip().startswith("PERCENT"):
                 cumul_balance_error = float(line.split()[3])
-                assert (
-                    abs(cumul_balance_error) < 0.00001
-                ), "Cumulative balance error = {} for {}, should equal 0.0".format(
-                    cumul_balance_error, mname
+                assert abs(cumul_balance_error) < 0.00001, (
+                    f"Cumulative balance error = {cumul_balance_error} for {mname}, "
+                    "should equal 0.0"
                 )
 
 
diff --git a/autotest/test_gwf_ifmod_idomain.py b/autotest/test_gwf_ifmod_idomain.py
index 18e3d9c6354..206986ce571 100644
--- a/autotest/test_gwf_ifmod_idomain.py
+++ b/autotest/test_gwf_ifmod_idomain.py
@@ -350,10 +350,8 @@ def check_output(idx, test):
 
         # compare heads
         maxdiff = np.amax(abs(heads - heads_2models))
-        assert maxdiff < 10 * hclose_check, "Max. head diff. {} should \
-                         be within solver tolerance (x10): {}".format(
-            maxdiff, 10 * hclose_check
-        )
+        assert maxdiff < 10 * hclose_check, f"Max. head diff. {maxdiff} should \
+                         be within solver tolerance (x10): {10 * hclose_check}"
 
     # check budget error from .lst file
     for mname in [mname_ref, mname_left, mname_right]:
@@ -361,10 +359,9 @@ def check_output(idx, test):
         for line in open(fpth):
             if line.lstrip().startswith("PERCENT"):
                 cumul_balance_error = float(line.split()[3])
-                assert (
-                    abs(cumul_balance_error) < 0.00001
-                ), "Cumulative balance error = {} for {}, should equal 0.0".format(
-                    cumul_balance_error, mname
+                assert abs(cumul_balance_error) < 0.00001, (
+                    f"Cumulative balance error = {cumul_balance_error} for {mname}, "
+                    "should equal 0.0"
                 )
 
 
diff --git a/autotest/test_gwf_ifmod_newton.py b/autotest/test_gwf_ifmod_newton.py
index 64b98e73872..bc6d46a9270 100644
--- a/autotest/test_gwf_ifmod_newton.py
+++ b/autotest/test_gwf_ifmod_newton.py
@@ -362,10 +362,8 @@ def check_output(idx, test):
 
         # compare heads
         maxdiff = np.amax(abs(heads - heads_2models))
-        assert maxdiff < 10 * hclose_check, "Max. head diff. {} should \
-                         be within solver tolerance (x10): {}".format(
-            maxdiff, 10 * hclose_check
-        )
+        assert maxdiff < 10 * hclose_check, f"Max. head diff. {maxdiff} should \
+                         be within solver tolerance (x10): {10 * hclose_check}"
 
     # check budget error from .lst file
     for mname in [mname_ref, mname_left, mname_right]:
@@ -373,10 +371,9 @@ def check_output(idx, test):
         for line in open(fpth):
             if line.lstrip().startswith("PERCENT"):
                 cumul_balance_error = float(line.split()[3])
-                assert (
-                    abs(cumul_balance_error) < 0.00001
-                ), "Cumulative balance error = {} for {}, should equal 0.0".format(
-                    cumul_balance_error, mname
+                assert abs(cumul_balance_error) < 0.00001, (
+                    f"Cumulative balance error = {cumul_balance_error} for {mname}, "
+                    "should equal 0.0"
                 )
 
 
diff --git a/autotest/test_gwf_ifmod_rewet.py b/autotest/test_gwf_ifmod_rewet.py
index d69c6e6aac5..f213a25bf37 100644
--- a/autotest/test_gwf_ifmod_rewet.py
+++ b/autotest/test_gwf_ifmod_rewet.py
@@ -383,10 +383,8 @@ def check_output(idx, test):
 
         # compare heads
         maxdiff = np.amax(abs(heads - heads_2models))
-        assert maxdiff < 10 * hclose_check, "Max. head diff. {} should \
-                         be within solver tolerance (x10): {}".format(
-            maxdiff, 10 * hclose_check
-        )
+        assert maxdiff < 10 * hclose_check, f"Max. head diff. {maxdiff} should \
+                         be within solver tolerance (x10): {10 * hclose_check}"
 
     # check budget error from .lst file
     for mname in [mname_ref, mname_left, mname_right]:
@@ -394,10 +392,9 @@ def check_output(idx, test):
         for line in open(fpth):
             if line.lstrip().startswith("PERCENT"):
                 cumul_balance_error = float(line.split()[3])
-                assert (
-                    abs(cumul_balance_error) < 0.00001
-                ), "Cumulative balance error = {} for {}, should equal 0.0".format(
-                    cumul_balance_error, mname
+                assert abs(cumul_balance_error) < 0.00001, (
+                    f"Cumulative balance error = {cumul_balance_error} for {mname}, "
+                    "should equal 0.0"
                 )
 
 
diff --git a/autotest/test_gwf_ifmod_xt3d01.py b/autotest/test_gwf_ifmod_xt3d01.py
index 08b0193b9f5..f9546b2492b 100644
--- a/autotest/test_gwf_ifmod_xt3d01.py
+++ b/autotest/test_gwf_ifmod_xt3d01.py
@@ -334,12 +334,10 @@ def exact(x):
             if h != 1e30:
                 diff = abs(h - exact(xc))
                 assert diff < 10 * hclose, (
-                    "head difference in parent model {}"
-                    " exceeds solver tolerance (x10) {}"
-                    " for row {} and col {}\n"
-                    "(should be {}, was {})".format(
-                        diff, 10 * hclose, irow + 1, icol + 1, exact(xc), h
-                    )
+                    f"head difference in parent model {diff}"
+                    f" exceeds solver tolerance (x10) {10 * hclose}"
+                    f" for row {irow + 1} and col {icol + 1}\n"
+                    f"(should be {exact(xc)}, was {h})"
                 )
 
     for irow in range(mg.nrow):
@@ -352,11 +350,9 @@ def exact(x):
             if qxb[0, irow, icol] != "--":
                 diff = abs(qxb[0, irow, icol] - qx_exact)
                 assert diff < 10 * hclose, (
-                    "Difference in spec. dis. for parent {}"
-                    " exceeds solver tolerance (x10) {}"
-                    " for row {} and col {}".format(
-                        diff, 10 * hclose, irow + 1, icol + 1
-                    )
+                    f"Difference in spec. dis. for parent {diff}"
+                    f" exceeds solver tolerance (x10) {10 * hclose}"
+                    f" for row {irow + 1} and col {icol + 1}"
                 )
 
     # and now the child
@@ -367,11 +363,9 @@ def exact(x):
             if h != 1e30:
                 diff = abs(h - exact(xc))
                 assert diff < 10 * hclose, (
-                    "Head difference in child model {}"
-                    " exceeds solver tolerance (x10) {}"
-                    " for row {} and col {}".format(
-                        diff, 10 * hclose, irow + 1, icol + 1
-                    )
+                    f"Head difference in child model {diff}"
+                    f" exceeds solver tolerance (x10) {10 * hclose}"
+                    f" for row {irow + 1} and col {icol + 1}"
                 )
 
     for irow in range(mg_c.nrow):
@@ -384,11 +378,9 @@ def exact(x):
             if qxb_c[0, irow, icol] != "--":
                 diff = abs(qxb_c[0, irow, icol] - qx_exact)
                 assert diff < 10 * hclose, (
-                    "Difference in spec. dis. for child {}"
-                    " exceeds solver tolerance (x10) {}"
-                    " for row {} and col {}".format(
-                        diff, 10 * hclose, irow + 1, icol + 1
-                    )
+                    f"Difference in spec. dis. for child {diff}"
+                    f" exceeds solver tolerance (x10) {10 * hclose}"
+                    f" for row {irow + 1} and col {icol + 1}"
                 )
 
     # todo: mflistbudget
@@ -398,10 +390,9 @@ def exact(x):
         for line in open(fpth):
             if line.lstrip().startswith("PERCENT"):
                 cumul_balance_error = float(line.split()[3])
-                assert (
-                    abs(cumul_balance_error) < 0.00001
-                ), "Cumulative balance error = {} for {}, should equal 0.0".format(
-                    cumul_balance_error, mname
+                assert abs(cumul_balance_error) < 0.00001, (
+                    f"Cumulative balance error = {cumul_balance_error} for {mname}, "
+                    "should equal 0.0"
                 )
 
     # Check on residual, which is stored in diagonal position of
diff --git a/autotest/test_gwf_ifmod_xt3d02.py b/autotest/test_gwf_ifmod_xt3d02.py
index d72aa9ef358..e9e9bb744ed 100644
--- a/autotest/test_gwf_ifmod_xt3d02.py
+++ b/autotest/test_gwf_ifmod_xt3d02.py
@@ -334,52 +334,38 @@ def check_output(idx, test):
 
     # compare heads
     maxdiff = np.amax(abs(heads - heads_2models))
-    assert maxdiff < 10 * hclose_check, "Max. head diff. {} should \
-                     be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. head diff. {maxdiff} should \
+                     be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_x left
     maxdiff = np.amax(abs(qxb[:, :, 0:5] - qxb_left))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (x) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (x) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_y left
     maxdiff = np.amax(abs(qyb[:, :, 0:5] - qyb_left))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (y) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (y) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_z left
     maxdiff = np.amax(abs(qzb[:, :, 0:5] - qzb_left))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (z) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (z) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_x right
     maxdiff = np.amax(abs(qxb[:, :, 5:] - qxb_right))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (x) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (x) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_y right
     maxdiff = np.amax(abs(qyb[:, :, 5:] - qyb_right))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (y) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (y) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_z right
     maxdiff = np.amax(abs(qzb[:, :, 5:] - qzb_right))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (z) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (z) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # check budget error from .lst file
     for mname in [mname_ref, mname_left, mname_right]:
@@ -387,10 +373,9 @@ def check_output(idx, test):
         for line in open(fpth):
             if line.lstrip().startswith("PERCENT"):
                 cumul_balance_error = float(line.split()[3])
-                assert (
-                    abs(cumul_balance_error) < 0.00001
-                ), "Cumulative balance error = {} for {}, should equal 0.0".format(
-                    cumul_balance_error, mname
+                assert abs(cumul_balance_error) < 0.00001, (
+                    f"Cumulative balance error = {cumul_balance_error} for {mname}, "
+                    "should equal 0.0"
                 )
 
 
diff --git a/autotest/test_gwf_ifmod_xt3d03.py b/autotest/test_gwf_ifmod_xt3d03.py
index 74f7ee1021f..d048efed2e6 100644
--- a/autotest/test_gwf_ifmod_xt3d03.py
+++ b/autotest/test_gwf_ifmod_xt3d03.py
@@ -370,10 +370,8 @@ def check_output(idx, test):
 
     # compare heads
     maxdiff = np.amax(abs(heads - heads_merged))
-    assert maxdiff < 10 * hclose_check, "Max. head diff. {} should \
-                     be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. head diff. {maxdiff} should \
+                     be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis-x
     qx_top = np.append(qx_tl[0], qx_tr[0], axis=1)
@@ -381,10 +379,8 @@ def check_output(idx, test):
     qx_merged = np.append(qx_top, qx_bot, axis=0)
 
     maxdiff = np.amax(abs(qx - qx_merged))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (x) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (x) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis-y
     qy_top = np.append(qy_tl[0], qy_tr[0], axis=1)
@@ -392,10 +388,8 @@ def check_output(idx, test):
     qy_merged = np.append(qy_top, qy_bot, axis=0)
 
     maxdiff = np.amax(abs(qy - qy_merged))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (y) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (y) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis-z
     qz_top = np.append(qz_tl[0], qz_tr[0], axis=1)
@@ -403,10 +397,8 @@ def check_output(idx, test):
     qz_merged = np.append(qz_top, qz_bot, axis=0)
 
     maxdiff = np.amax(abs(qz - qz_merged))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (z) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (z) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # check budget error from .lst file
     for mname in ["ref", "tl", "tr", "bl", "br"]:
@@ -414,10 +406,9 @@ def check_output(idx, test):
         for line in open(fpth):
             if line.lstrip().startswith("PERCENT"):
                 cumul_balance_error = float(line.split()[3])
-                assert (
-                    abs(cumul_balance_error) < 0.00001
-                ), "Cumulative balance error = {} for {}, should equal 0.0".format(
-                    cumul_balance_error, mname
+                assert abs(cumul_balance_error) < 0.00001, (
+                    f"Cumulative balance error = {cumul_balance_error} for {mname}, "
+                    "should equal 0.0"
                 )
 
 
diff --git a/autotest/test_gwf_maw10.py b/autotest/test_gwf_maw10.py
index 92ea7e6a5f2..4e41e234d37 100644
--- a/autotest/test_gwf_maw10.py
+++ b/autotest/test_gwf_maw10.py
@@ -245,7 +245,7 @@ def check_output(idx, test):
         msgtime = (
             "There should be a matching time in the maw.obs.csv file for each "
             "time in the maw.reduce.csv file, but no match was found for "
-            "time = {} in the maw.obs.csv file".format(timevalmfr)
+            f"time = {timevalmfr} in the maw.obs.csv file"
         )
         blnFoundTimeMatch = False
         for rowobs in tcobs:
@@ -254,9 +254,9 @@ def check_output(idx, test):
                 blnFoundTimeMatch = True
                 actvalobs = rowobs["M1RATE"]
                 msgval = (
-                    "The maw.obs.csv file rate-actual value of {} should have "
-                    "matched the maw.reduce.csv file rate-actual value of {} "
-                    "at time {}".format(actvalobs, actvalmfr, timevalobs)
+                    f"The maw.obs.csv file rate-actual value of {actvalobs} "
+                    "should have matched the maw.reduce.csv file rate-actual "
+                    f"value of {actvalmfr} at time {timevalobs}"
                 )
                 break
         assert blnFoundTimeMatch, msgtime
diff --git a/autotest/test_gwf_maw11.py b/autotest/test_gwf_maw11.py
index 0d48f12d615..f96bfaea496 100644
--- a/autotest/test_gwf_maw11.py
+++ b/autotest/test_gwf_maw11.py
@@ -188,7 +188,7 @@ def build_models(idx, test):
         connectiondata=connectiondata,
         perioddata=sfr_perioddata,
         pname="SFR-1",
-        filename="{}.sfr".format(name),
+        filename=f"{name}.sfr",
     )
 
     packages = [("MAW-1",), ("SFR-1",)]
@@ -202,7 +202,7 @@ def build_models(idx, test):
         packages=packages,
         perioddata=mvr_perioddata,
         pname="MVR-1",
-        filename="{}.mvr".format(name),
+        filename=f"{name}.mvr",
     )
 
     # output control
@@ -235,9 +235,9 @@ def check_output(idx, test):
     for itm in bud_names:
         nm_lst.append(str(itm.strip(), "utf-8"))
 
-    assert "CONSTANT-TO-MVR" in nm_lst, (
-        "Expected budget term not in MAW " "binary output file."
-    )
+    assert (
+        "CONSTANT-TO-MVR" in nm_lst
+    ), "Expected budget term not in MAW binary output file."
 
     rtm = mawobj.get_data(text="RATE-TO-MVR")
     ctm = mawobj.get_data(text="CONSTANT-TO-MVR")
diff --git a/autotest/test_gwf_mf6io_app2_examples_distypes.py b/autotest/test_gwf_mf6io_app2_examples_distypes.py
index f5a65ad6930..7c70effe231 100644
--- a/autotest/test_gwf_mf6io_app2_examples_distypes.py
+++ b/autotest/test_gwf_mf6io_app2_examples_distypes.py
@@ -266,7 +266,7 @@ def build_mf6(idx, ws, gridgen):
     elif "disu" in str(ws):
         dis_type = "disu"
     else:
-        raise ValueError(f"Invalid discretization type in {str(ws)}")
+        raise ValueError(f"Invalid discretization type in {ws!s}")
 
     if "disu" in str(ws):
         list_recharge = True
diff --git a/autotest/test_gwf_npf_tvk04.py b/autotest/test_gwf_npf_tvk04.py
index 4a01667e521..46b1e2d5489 100644
--- a/autotest/test_gwf_npf_tvk04.py
+++ b/autotest/test_gwf_npf_tvk04.py
@@ -171,8 +171,8 @@ def check_output(idx, test):
         sp_x.append(bud)
 
     # comment when done testing
-    print(f"Total outflow in stress period 1 is {str(sp_x[0][8])}")
-    print(f"Total outflow in stress period 2 after increasing K33 is {str(sp_x[1][8])}")
+    print(f"Total outflow in stress period 1 is {sp_x[0][8]!s}")
+    print(f"Total outflow in stress period 2 after increasing K33 is {sp_x[1][8]!s}")
     errmsg = (
         "Expect higher flow rate in period 2 compared to period 1, "
         "but found equal or higher flow rate in period 1"
diff --git a/autotest/test_gwf_npf_tvk05.py b/autotest/test_gwf_npf_tvk05.py
index 19097088e19..28723d5819d 100644
--- a/autotest/test_gwf_npf_tvk05.py
+++ b/autotest/test_gwf_npf_tvk05.py
@@ -172,7 +172,7 @@ def check_output(idx, test):
         sp_x.append(bud)
 
     # comment when done testing
-    print(f"Total outflow in stress period 1 is {str(sp_x[0][8])}")
+    print(f"Total outflow in stress period 1 is {sp_x[0][8]!s}")
     print(
         "Total outflow in stress period 2 after increasing K33 "
         "should have no effect on final solution"
diff --git a/autotest/test_gwf_returncodes.py b/autotest/test_gwf_returncodes.py
index 20abe43f51f..9b60a476a76 100644
--- a/autotest/test_gwf_returncodes.py
+++ b/autotest/test_gwf_returncodes.py
@@ -158,7 +158,7 @@ def converge_fail_continue(dir, exe):
     msg = (
         "The run should have been successful even though it failed, because"
         " the continue flag was set.  But a non-zero error code was "
-        "found: {}".format(returncode)
+        f"found: {returncode}"
     )
     assert returncode == 0, msg
 
diff --git a/autotest/test_gwf_uzf_gwet_opt.py b/autotest/test_gwf_uzf_gwet_opt.py
index 647373fe4c3..82fc08e76db 100644
--- a/autotest/test_gwf_uzf_gwet_opt.py
+++ b/autotest/test_gwf_uzf_gwet_opt.py
@@ -374,7 +374,7 @@ def build_gwf_model(sim, name, linear_et_flg=True):
         outer_dvclose=1e-5,
         inner_dvclose=1e-6,
         linear_acceleration="BICGSTAB",
-        filename="{}.ims".format(name),
+        filename=f"{name}.ims",
     )
     sim.register_ims_package(imsgwf, [gwf.name])
 
diff --git a/autotest/test_gwt_adv01_gwtgwt.py b/autotest/test_gwt_adv01_gwtgwt.py
index d00f5dd2108..b8e488f18d7 100644
--- a/autotest/test_gwt_adv01_gwtgwt.py
+++ b/autotest/test_gwt_adv01_gwtgwt.py
@@ -680,10 +680,9 @@ def check_output(idx, test):
         for line in open(fpth):
             if line.lstrip().startswith("PERCENT"):
                 cumul_balance_error = float(line.split()[3])
-                assert (
-                    abs(cumul_balance_error) < 0.00001
-                ), "Cumulative balance error = {} for {}, should equal 0.0".format(
-                    cumul_balance_error, mname
+                assert abs(cumul_balance_error) < 0.00001, (
+                    f"Cumulative balance error = {cumul_balance_error} for {mname}, "
+                    "should equal 0.0"
                 )
 
         # get grid data (from GWF)
diff --git a/autotest/test_gwt_henry_gwtgwt.py b/autotest/test_gwt_henry_gwtgwt.py
index 3d747951dc9..91488b81e79 100644
--- a/autotest/test_gwt_henry_gwtgwt.py
+++ b/autotest/test_gwt_henry_gwtgwt.py
@@ -357,8 +357,8 @@ def check_output(idx, test):
 
     # compare heads
     maxdiff = np.amax(abs(heads - heads_gwfgwf))
-    assert maxdiff < 10 * hclose, "Max. head diff. {} should \
-        be within solver tolerance (x10): {}".format(maxdiff, 10 * hclose)
+    assert maxdiff < 10 * hclose, f"Max. head diff. {maxdiff} should \
+        be within solver tolerance (x10): {10 * hclose}"
 
     fpth = os.path.join(test.workspace, "gwt_ref.ucn")
     try:
@@ -385,8 +385,8 @@ def check_output(idx, test):
     conc_gwtgwt = np.append(conc_left, conc_right, axis=2)
 
     maxdiff = np.amax(abs(conc_gwtgwt - conc_ref))
-    assert maxdiff < conc_tol, "Max. concentration diff. {} should \
-        be within solver tolerance (x10): {}".format(maxdiff, conc_tol)
+    assert maxdiff < conc_tol, f"Max. concentration diff. {maxdiff} should \
+        be within solver tolerance (x10): {conc_tol}"
 
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))
diff --git a/autotest/test_gwt_moc3d01_zod.py b/autotest/test_gwt_moc3d01_zod.py
index 43af4756229..fe705666145 100644
--- a/autotest/test_gwt_moc3d01_zod.py
+++ b/autotest/test_gwt_moc3d01_zod.py
@@ -400,8 +400,8 @@ def check_output(idx, test):
         if conc[i] / delt > decay_rate:
             qknown = -decay_rate * vcell * porosity
             errmsg = (
-                "Decay rate in budget file for cell {} should be "
-                "{} but found {} instead.".format(i, qdecay_budfile[i], qknown)
+                f"Decay rate in budget file for cell {i} should be "
+                f"{qdecay_budfile[i]} but found {qknown} instead."
             )
             assert np.allclose(qdecay_budfile[i], qknown), errmsg
         # print(i, qdecay_budfile[i], conc[i])
@@ -427,7 +427,7 @@ def check_output(idx, test):
             "Mass transfer rates from the gwt budget file do not "
             "compare with mass transfer rates calculated from "
             "simulated mobile and immobile domain concentrations\n"
-            "{} /= {}".format(qim_budfile, qim_calculated)
+            f"{qim_budfile} /= {qim_calculated}"
         )
         np.allclose(qim_budfile, qim_calculated), errmsg
 
@@ -556,7 +556,7 @@ def check_output(idx, test):
     tsres = tsreslist[idx]
     errmsg = (
         "Simulated concentrations do not match with known solution.\n"
-        "{} /= {}".format(tssim, tsres)
+        f"{tssim} /= {tsres}"
     )
     if tsres is not None:
         assert np.allclose(tsres, tssim), errmsg
diff --git a/autotest/test_gwt_mst06_noadv.py b/autotest/test_gwt_mst06_noadv.py
index 07f88395ba6..b583d84f7c6 100644
--- a/autotest/test_gwt_mst06_noadv.py
+++ b/autotest/test_gwt_mst06_noadv.py
@@ -122,10 +122,10 @@ def check_output(idx, test):
     # The answer
     # print(conc[:, 1])
     cres = np.array([7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0, 0.0, 0.0])
-    msg = "simulated concentrations do not match with known solution. {} {}".format(
-        conc[:, 1], cres
+    assert np.allclose(cres, conc[:, 1]), (
+        "simulated concentrations do not match with known solution. "
+        f"{conc[:, 1]} {cres}"
     )
-    assert np.allclose(cres, conc[:, 1]), msg
 
     # Check budget file
     fpth = os.path.join(test.workspace, f"{gwtname}.bud")
@@ -136,10 +136,10 @@ def check_output(idx, test):
     decay_list = bobj.get_data(text="DECAY-AQUEOUS")
     decay_rate = [dr[0] for dr in decay_list]
     decay_rate_answer = [-8.4, -8.4, -8.4, -8.4, -8.4, -8.4, -8.4, -8.4, 0.0, 0.0]
-    msg = "simulated decay rates do not match with known solution. {} {}".format(
-        decay_rate, decay_rate_answer
+    assert np.allclose(decay_rate, decay_rate_answer), (
+        "simulated decay rates do not match with known solution. "
+        f"{decay_rate} {decay_rate_answer}"
     )
-    assert np.allclose(decay_rate, decay_rate_answer), msg
 
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))
diff --git a/autotest/test_gwt_mt3dms_p01.py b/autotest/test_gwt_mt3dms_p01.py
index f8875b73a6b..99e9429f248 100644
--- a/autotest/test_gwt_mt3dms_p01.py
+++ b/autotest/test_gwt_mt3dms_p01.py
@@ -20,7 +20,6 @@
 
 import os
 from pathlib import Path
-from typing import Tuple
 
 import flopy
 import numpy as np
@@ -461,7 +460,7 @@ def p01mf6(
     return sim, conc
 
 
-def get_binaries(targets) -> Tuple[Path, Path, Path]:
+def get_binaries(targets) -> tuple[Path, Path, Path]:
     return (
         targets["mf6"],
         try_get_target(targets, "mf2005s"),
diff --git a/autotest/test_gwt_sft_inactive01.py b/autotest/test_gwt_sft_inactive01.py
index b10ad2ac28e..4743f0a1265 100644
--- a/autotest/test_gwt_sft_inactive01.py
+++ b/autotest/test_gwt_sft_inactive01.py
@@ -388,7 +388,7 @@ def build_models(idx, test):
         outer_dvclose=1e-5,
         inner_dvclose=1e-6,
         linear_acceleration="BICGSTAB",
-        filename="{}.ims".format(name),
+        filename=f"{name}.ims",
     )
     sim.register_ims_package(imsgwf, [gwf.name])
 
@@ -441,14 +441,14 @@ def build_models(idx, test):
         auxiliary="CONCENTRATION",
         stress_period_data=chd_spd_left,
         pname="WEL-left",
-        filename="{}.wel-left".format(name),
+        filename=f"{name}.wel-left",
     )
     flopy.mf6.ModflowGwfchd(
         gwf,
         auxiliary="CONCENTRATION",
         stress_period_data=chd_spd_right,
         pname="WEL-right",
-        filename="{}.wel-right".format(name),
+        filename=f"{name}.wel-right",
     )
 
     # SFR data
@@ -560,7 +560,7 @@ def build_models(idx, test):
         outer_dvclose=1e-5,
         inner_dvclose=1e-6,
         linear_acceleration="BICGSTAB",
-        filename="{}.ims".format(gwtname),
+        filename=f"{gwtname}.ims",
     )
     sim.register_ims_package(imsgwt, [gwt.name])
 
diff --git a/autotest/test_gwtgwt_oldexg.py b/autotest/test_gwtgwt_oldexg.py
index b71399982c5..7ed5ec5e266 100644
--- a/autotest/test_gwtgwt_oldexg.py
+++ b/autotest/test_gwtgwt_oldexg.py
@@ -605,52 +605,38 @@ def compare_gwf_to_ref(test):
 
     # compare heads
     maxdiff = np.amax(abs(heads - heads_2models))
-    assert maxdiff < 10 * hclose_check, "Max. head diff. {} should \
-                     be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. head diff. {maxdiff} should \
+                     be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_x left
     maxdiff = np.amax(abs(qxb[:, :, 0:5] - qxb_left))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (x) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (x) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_y left
     maxdiff = np.amax(abs(qyb[:, :, 0:5] - qyb_left))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (y) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (y) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_z left
     maxdiff = np.amax(abs(qzb[:, :, 0:5] - qzb_left))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (z) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (z) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_x right
     maxdiff = np.amax(abs(qxb[:, :, 5:] - qxb_right))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (x) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (x) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_y right
     maxdiff = np.amax(abs(qyb[:, :, 5:] - qyb_right))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (y) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (y) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # compare spdis_z right
     maxdiff = np.amax(abs(qzb[:, :, 5:] - qzb_right))
-    assert maxdiff < 10 * hclose_check, "Max. diff. in spec. discharge (z) {} \
-                     should be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. diff. in spec. discharge (z) {maxdiff} \
+                     should be within solver tolerance (x10): {10 * hclose_check}"
 
     # check budget error from .lst file
     for mname in [mname_ref, mname_left, mname_right]:
@@ -658,10 +644,9 @@ def compare_gwf_to_ref(test):
         for line in open(fpth):
             if line.lstrip().startswith("PERCENT"):
                 cumul_balance_error = float(line.split()[3])
-                assert (
-                    abs(cumul_balance_error) < 0.00001
-                ), "Cumulative balance error = {} for {}, should equal 0.0".format(
-                    cumul_balance_error, mname
+                assert abs(cumul_balance_error) < 0.00001, (
+                    f"Cumulative balance error = {cumul_balance_error} for {mname}, "
+                    "should equal 0.0"
                 )
 
     # check flowja residual
@@ -704,10 +689,8 @@ def compare_gwt_to_ref(test):
 
     # compare concentrations
     maxdiff = np.amax(abs(conc - conc_2models))
-    assert maxdiff < 10 * hclose_check, "Max. concentration diff. {} should \
-                     be within solver tolerance (x10): {}".format(
-        maxdiff, 10 * hclose_check
-    )
+    assert maxdiff < 10 * hclose_check, f"Max. concentration diff. {maxdiff} should \
+                     be within solver tolerance (x10): {10 * hclose_check}"
 
     # check budget error from .lst file
     for mname in [mname_gwtref, mname_gwtleft, mname_gwtright]:
@@ -715,10 +698,9 @@ def compare_gwt_to_ref(test):
         for line in open(fpth):
             if line.lstrip().startswith("PERCENT"):
                 cumul_balance_error = float(line.split()[3])
-                assert (
-                    abs(cumul_balance_error) < 0.00001
-                ), "Cumulative balance error = {} for {}, should equal 0.0".format(
-                    cumul_balance_error, mname
+                assert abs(cumul_balance_error) < 0.00001, (
+                    f"Cumulative balance error = {cumul_balance_error} for {mname}, "
+                    "should equal 0.0"
                 )
 
     # check flowja residual for transport mass flows
diff --git a/autotest/test_olf_dis.py b/autotest/test_olf_dis.py
index 0c9e9994bb5..c27eaa9dd66 100644
--- a/autotest/test_olf_dis.py
+++ b/autotest/test_olf_dis.py
@@ -229,8 +229,8 @@ def check_grb_dis2d(fpth):
     assert grb.xorigin == xorigin, "xorigin in grb file is not correct"
     assert grb.yorigin == yorigin, "yorigin in grb file is not correct"
     assert grb.angrot == angrot, "angrot in grb file is not correct"
-    assert np.allclose(grb.delr, dx * np.ones((ncol))), "grb delr not correct"
-    assert np.allclose(grb.delc, dx * np.ones((nrow))), "grb delc not correct"
+    assert np.allclose(grb.delr, dx * np.ones(ncol)), "grb delr not correct"
+    assert np.allclose(grb.delc, dx * np.ones(nrow)), "grb delc not correct"
     assert np.allclose(
         grb.bot.reshape((nrow, ncol)), np.zeros((nrow, ncol))
     ), "grb botm not correct"
diff --git a/autotest/test_prt_track_events.py b/autotest/test_prt_track_events.py
index 8c28f0f18df..4de9bd5f8f2 100644
--- a/autotest/test_prt_track_events.py
+++ b/autotest/test_prt_track_events.py
@@ -23,7 +23,6 @@
 """
 
 from pathlib import Path
-from typing import List
 
 import flopy
 import matplotlib.cm as cm
@@ -147,7 +146,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         for grp in ["a", "b"]
     ]
 
-    def get_oc() -> List[str]:
+    def get_oc() -> list[str]:
         prt_track_file = f"{prt_name}.trk"
         prt_track_csv_file = f"{prt_name}.trk.csv"
         if "all" in name:
diff --git a/distribution/benchmark.py b/distribution/benchmark.py
index 203898b246a..7c335b94e2b 100644
--- a/distribution/benchmark.py
+++ b/distribution/benchmark.py
@@ -7,7 +7,6 @@
 from multiprocessing import Pool
 from os import PathLike
 from pathlib import Path
-from typing import List, Tuple
 
 import flopy
 import pytest
@@ -34,7 +33,7 @@
 )
 
 
-def download_previous_version(output_path: PathLike) -> Tuple[str, Path]:
+def download_previous_version(output_path: PathLike) -> tuple[str, Path]:
     output_path = Path(output_path).expanduser().absolute()
     version = get_latest_version(GITHUB_REPO)
     distname = f"mf{version}_{OSTAG}"
@@ -258,7 +257,7 @@ def write_results(
     output_path: PathLike,
     current_total,
     previous_total,
-    lines: List[str],
+    lines: list[str],
 ):
     current_exe = Path(current_exe)
     previous_exe = Path(previous_exe)
@@ -313,7 +312,7 @@ def run_benchmarks(
     previous_bin_path: PathLike,
     examples_path: PathLike,
     output_path: PathLike,
-    excluded: List[str] = [],
+    excluded: list[str] = [],
 ):
     """Benchmark current development version against previous release
     with example models.
diff --git a/distribution/build_dist.py b/distribution/build_dist.py
index b103896509b..2f46ef86441 100644
--- a/distribution/build_dist.py
+++ b/distribution/build_dist.py
@@ -6,7 +6,7 @@
 from pathlib import Path
 from pprint import pprint
 from shutil import copy, copyfile, copytree, ignore_patterns, rmtree
-from typing import List, Optional
+from typing import Optional
 
 import pytest
 from build_docs import build_documentation
@@ -110,7 +110,7 @@ def setup_examples(
     bin_path: PathLike,
     examples_path: PathLike,
     force: bool = False,
-    models: Optional[List[str]] = None,
+    models: Optional[list[str]] = None,
 ):
     examples_path = Path(examples_path).expanduser().absolute()
 
@@ -285,7 +285,7 @@ def build_distribution(
     output_path: PathLike,
     full: bool = False,
     force: bool = False,
-    models: Optional[List[str]] = None,
+    models: Optional[list[str]] = None,
 ):
     print(f"Building {'full' if full else 'minimal'} distribution")
 
diff --git a/distribution/build_docs.py b/distribution/build_docs.py
index 4ea33a1c459..5ae0fa4fec1 100644
--- a/distribution/build_docs.py
+++ b/distribution/build_docs.py
@@ -9,7 +9,7 @@
 from pathlib import Path
 from pprint import pprint
 from tempfile import TemporaryDirectory
-from typing import List, Optional
+from typing import Optional
 from urllib.error import HTTPError
 from warnings import warn
 
@@ -210,7 +210,7 @@ def test_build_deprecations_tex():
     build_deprecations_tex(force=True)
 
 
-def build_mf6io_tex(models: Optional[List[str]] = None, force: bool = False):
+def build_mf6io_tex(models: Optional[list[str]] = None, force: bool = False):
     """Build LaTeX files for the MF6IO guide from DFN files."""
 
     if models is None:
@@ -326,7 +326,7 @@ def build_usage_example_tex(
 
 
 def build_pdfs(
-    tex_paths: List[PathLike],
+    tex_paths: list[PathLike],
     output_path: PathLike,
     passes: int = 3,
     force: bool = False,
@@ -401,7 +401,7 @@ def build_documentation(
     output_path: PathLike,
     force: bool = False,
     full: bool = False,
-    models: Optional[List[str]] = None,
+    models: Optional[list[str]] = None,
     repo_owner: str = "MODFLOW-USGS",
 ):
     """Build documentation for a MODFLOW 6 distribution."""
diff --git a/distribution/conftest.py b/distribution/conftest.py
index 3165c8c04f5..024b6cedd6d 100644
--- a/distribution/conftest.py
+++ b/distribution/conftest.py
@@ -5,7 +5,7 @@
 PROJ_ROOT_PATH = Path(__file__).resolve().parent.parent
 DIST_PATH = (
     PROJ_ROOT_PATH.parent
-    / f"mf{str(Version((PROJ_ROOT_PATH / 'version.txt').read_text().strip()))}"
+    / f"mf{Version((PROJ_ROOT_PATH / 'version.txt').read_text().strip())!s}"
 )
 
 
diff --git a/distribution/utils.py b/distribution/utils.py
index bc79715b1fd..ef0bbc1b6b1 100644
--- a/distribution/utils.py
+++ b/distribution/utils.py
@@ -1,10 +1,11 @@
 import re
 import subprocess
 import sys
+from collections.abc import Iterator
 from datetime import datetime
 from pathlib import Path
 from pprint import pformat
-from typing import Iterator, List, Optional
+from typing import Optional
 
 _project_root_path = Path(__file__).resolve().parent.parent
 
@@ -20,8 +21,8 @@ def get_modified_time(path: Path) -> float:
 def glob(
     path: Path,
     pattern: str,
-    included: Optional[List[str]],
-    excluded: Optional[List[str]],
+    included: Optional[list[str]],
+    excluded: Optional[list[str]],
 ) -> Iterator[Path]:
     def is_included(p):
         if included is None:
diff --git a/doc/SuppTechInfo/python/STO-SpecificStorage.ipynb b/doc/SuppTechInfo/python/STO-SpecificStorage.ipynb
index d58a63c90d5..86134cb4337 100644
--- a/doc/SuppTechInfo/python/STO-SpecificStorage.ipynb
+++ b/doc/SuppTechInfo/python/STO-SpecificStorage.ipynb
@@ -124,7 +124,7 @@
     "    (0.0, 0.1, 0.2, 0.3, 0.4, 0.5), (\"red\", \"orange\", \"green\", \"cyan\", \"blue\", \"black\")\n",
     "):\n",
     "    sf1 = 1.0 - sf0\n",
-    "    label = r\"S$_F^t$ = {:.2f} S$_F^{{told}}$ = {:.2f}\".format(sf1, sf0)\n",
+    "    label = rf\"S$_F^t$ = {sf1:.2f} S$_F^{{told}}$ = {sf0:.2f}\"\n",
     "    ax.plot(\n",
     "        top,\n",
     "        calc_err(top, ss=1, sf0=sf0, sf1=sf1),\n",
@@ -159,7 +159,7 @@
     "    (0.0, 0.1, 0.2, 0.3, 0.4, 0.5), (\"red\", \"orange\", \"green\", \"cyan\", \"blue\", \"black\")\n",
     "):\n",
     "    sf1 = 1.0 - sf0\n",
-    "    label = r\"S$_F^t$ = {:.2f} S$_F^{{told}}$ = {:.2f}\".format(sf1, sf0)\n",
+    "    label = rf\"S$_F^t$ = {sf1:.2f} S$_F^{{told}}$ = {sf0:.2f}\"\n",
     "    ax.plot(\n",
     "        top,\n",
     "        err_function(top, ss=1, sf0=sf0, sf1=sf1),\n",
@@ -191,7 +191,7 @@
     "    (0.0, 0.1, 0.2, 0.3, 0.4, 0.5), (\"red\", \"orange\", \"green\", \"cyan\", \"blue\", \"black\")\n",
     "):\n",
     "    sf1 = 1.0 - sf0\n",
-    "    label = r\"S$_F^t$ = {:.2f} S$_F^{{told}}$ = {:.2f}\".format(sf1, sf0)\n",
+    "    label = rf\"S$_F^t$ = {sf1:.2f} S$_F^{{told}}$ = {sf0:.2f}\"\n",
     "    ax.plot(\n",
     "        top,\n",
     "        err_function_simp(top, ss=1e-5, sf0=sf0, sf1=sf1),\n",
diff --git a/doc/mf6io/mf6ivar/deprecations.py b/doc/mf6io/mf6ivar/deprecations.py
index 3a8d89a42be..353d441ea37 100644
--- a/doc/mf6io/mf6ivar/deprecations.py
+++ b/doc/mf6io/mf6ivar/deprecations.py
@@ -1,5 +1,5 @@
 from pathlib import Path
-from typing import List, Optional, Tuple
+from typing import Optional
 
 from packaging.version import Version
 
@@ -9,7 +9,7 @@
 
 def get_deprecations(
     dfndir,
-) -> List[Tuple[Path, str, Version, Optional[Version]]]:
+) -> list[tuple[Path, str, Version, Optional[Version]]]:
     dfns = Path(dfndir).rglob("*.dfn")
     deps = {}
     for dfn in dfns:
diff --git a/doc/mf6io/mf6ivar/fortran_parser.py b/doc/mf6io/mf6ivar/fortran_parser.py
index 7a1797e0be2..0fbf4c6964e 100644
--- a/doc/mf6io/mf6ivar/fortran_parser.py
+++ b/doc/mf6io/mf6ivar/fortran_parser.py
@@ -5,7 +5,7 @@ def get_next_line(line_list):
     line = ""
     while len(line_list) > 0:
         line = line_list.pop(0)
-        line = r"{}".format(line)
+        line = rf"{line}"
         line = line.strip()
         if "!" in line:
             idx = line.index("!")
diff --git a/doc/mf6io/mf6ivar/mf6ivar.py b/doc/mf6io/mf6ivar/mf6ivar.py
index fae2434499a..34269d087df 100644
--- a/doc/mf6io/mf6ivar/mf6ivar.py
+++ b/doc/mf6io/mf6ivar/mf6ivar.py
@@ -339,7 +339,7 @@ def write_desc(vardict, block, blk_var_list, varexcludeprefix=None):
     for name, b in vardict.keys():
         v = vardict[(name, b)]
         if v["block"] == block:
-            if "block_variable" in v and v["block_variable"]:
+            if v.get("block_variable"):
                 optional = "optional" in v and v["optional"] == "true"
                 blk_var_list.append((v["name"], optional))
             addv = True
@@ -401,7 +401,7 @@ def write_desc_md(vardict, block, blk_var_list, varexcludeprefix=None):
     for name, b in vardict.keys():
         v = vardict[(name, b)]
         if v["block"] == block:
-            if "block_variable" in v and v["block_variable"]:
+            if v.get("block_variable"):
                 optional = "optional" in v and v["optional"] == "true"
                 blk_var_list.append((v["name"], optional))
             addv = True
@@ -652,10 +652,10 @@ def write_appendix(blocks):
                 and "time" in blockname.lower()
             ):
                 oc = "no"
-            s = "{} & {} & {} & {} \\\\ \n".format(
-                component.upper(), ftype.upper(), blockname.upper(), oc
+            f.write(
+                f"{component.upper()} & {ftype.upper()} & {blockname.upper()} & {oc} "
+                "\\\\ \n"
             )
-            f.write(s)
             lastftype = ftype
 
         f.write("\n\n\\hline\n\\end{longtable}\n\\label{table:blocks}\n\\normalsize\n")
diff --git a/ruff.toml b/ruff.toml
index c0a4134f39e..59979055ba5 100644
--- a/ruff.toml
+++ b/ruff.toml
@@ -7,11 +7,16 @@ select = [
     "E",    # pycodestyle error
     "F",    # Pyflakes
     "I001", # isort - unsorted-imports
+    "RUF",  # Ruff-specific rules
+    "UP",   # Pyupgrade
 ]
 ignore = [
     "E722", # do not use bare `except`
     "E741", # ambiguous variable name
     "F841", # local variable assigned but never used
+    "RUF005", # collection literal concatenation
+    "RUF012", # mutable class default
+    "UP015", # redundant open modes
 ]
 
 [lint.per-file-ignores]
diff --git a/utils/idmloader/scripts/dfn2f90.py b/utils/idmloader/scripts/dfn2f90.py
index 6478ca06b20..eb7b430f384 100644
--- a/utils/idmloader/scripts/dfn2f90.py
+++ b/utils/idmloader/scripts/dfn2f90.py
@@ -3,6 +3,7 @@
 import textwrap
 from pathlib import Path
 from pprint import pprint
+from typing import Optional
 
 import yaml
 
@@ -18,7 +19,7 @@
 class Dfn2F90:
     """generate idm f90 file from dfn file"""
 
-    def __init__(self, dfnfspec: str = None, verbose: bool = False):
+    def __init__(self, dfnfspec: Optional[str] = None, verbose: bool = False):
         """Dfn290 init"""
 
         self._dfnfspec = dfnfspec
@@ -555,8 +556,8 @@ class IdmDfnSelector:
 
     def __init__(
         self,
-        dfn_d: dict = None,
-        varnames: list = None,
+        dfn_d: Optional[dict] = None,
+        varnames: Optional[list] = None,
     ):
         """IdmDfnSelector init"""