diff --git a/.doc/conf.py b/.doc/conf.py
index b2149eda..150a857e 100644
--- a/.doc/conf.py
+++ b/.doc/conf.py
@@ -39,7 +39,7 @@
     lines += (len(lines) - 1) * "-" + "\n\n"
     lines += (
         "The Jupyter Notebooks used to create the input files and figures for \n"
-        + "each of the MODFLOW 6 `examples <examples.html>`_.\n\n"
+        "each of the MODFLOW 6 `examples <examples.html>`_.\n\n"
     )
     f.write(lines)
 
diff --git a/etc/ci_create_examples_rst.py b/etc/ci_create_examples_rst.py
index 34291645..4cad9c4e 100644
--- a/etc/ci_create_examples_rst.py
+++ b/etc/ci_create_examples_rst.py
@@ -112,7 +112,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"))
 
     # read restructured text file as a string
     print(f"reading...'{dst}'")
@@ -125,10 +126,10 @@
     eq_labels = []
     for v in ex_regex.findall(lines):
         tag = f"\\label{{{v}}}"
-        label = "   :label: {}".format(v.replace(":", "-"))
+        label = f"   :label: {v.replace(':', '-')}"
         eq_labels.append(label)
-        tag = "`[{0}] <#{0}>`__".format(v)
-        replace_eq_labels[tag] = ":eq:`{}`".format(v.replace(":", "-"))
+        tag = f"`[{v}] <#{v}>`__"
+        replace_eq_labels[tag] = f":eq:`{v.replace(':', '-')}`"
 
     # find figure references in lines
     ex_regex = re.compile("\\`(.*?)\\`__")
@@ -168,7 +169,7 @@
                 if key in line:
                     line = line.replace(key, value)
 
-            tag = " "
+            tag = "\xa0"  # non-breaking space
             if tag in line:
                 line = line.replace(tag, " ")
 
@@ -273,9 +274,9 @@
         line += "----------------\n\n"
         line += (
             "The Jupyter notebook used to create the MODFLOW 6 input files\n"
-            + "for this example and post-process the results is:\n\n"
+            "for this example and post-process the results is:\n\n"
         )
-        line += "* `{0} <../_notebooks/{0}.html>`_\n".format(ex_root)
+        line += f"* `{ex_root} <../_notebooks/{ex_root}.html>`_\n"
         line += "\n"
 
         # Check to see if there is a gif with the same name as the example name
diff --git a/ruff.toml b/ruff.toml
index bf6c4dbf..e2c19576 100644
--- a/ruff.toml
+++ b/ruff.toml
@@ -7,10 +7,15 @@ select = [
     "E",    # pycodestyle error
     "F",    # Pyflakes
     "I001", # isort - unsorted-imports
+    # "ISC",  # implicitly concatenated strings
+    "RUF",  # Ruff-specific rules
+    "UP",   # Pyupgrade
 ]
 ignore = [
     "E501", # line too long TODO FIXME
     "E722", # do not use bare `except`
     "E741", # ambiguous variable name
     "F841", # local variable assigned but never used
+    "RUF005", # collection literal concatenation
+    "UP015", # redundant open modes
 ]
diff --git a/scripts/ex-gwe-ates.py b/scripts/ex-gwe-ates.py
index 90c35b2c..ea193618 100644
--- a/scripts/ex-gwe-ates.py
+++ b/scripts/ex-gwe-ates.py
@@ -564,7 +564,7 @@ def build_model(sim_name, verts, cell2d, top, botm):
         vertices=verts,
         cell2d=cell2d,
         pname="DISV",
-        filename="{}.disv".format(gwfname),
+        filename=f"{gwfname}.disv",
     )
 
     # Instantiating MODFLOW 6 node property flow package
@@ -582,7 +582,7 @@ def build_model(sim_name, verts, cell2d, top, botm):
         save_specific_discharge=True,
         save_saturation=True,
         pname="NPF",
-        filename="{}.npf".format(gwfname),
+        filename=f"{gwfname}.npf",
     )
 
     # Instantiating MODFLOW 6 initial conditions package
@@ -596,7 +596,7 @@ def build_model(sim_name, verts, cell2d, top, botm):
         sy=sy,
         transient={0: True},
         pname="STO",
-        filename="{}.sto".format(gwfname),
+        filename=f"{gwfname}.sto",
     )
 
     # Instantiate WEL package
@@ -622,12 +622,12 @@ def build_model(sim_name, verts, cell2d, top, botm):
         save_flows=True,
         stress_period_data=wel_spd,
         pname="WEL",
-        filename="{}.wel".format(gwfname),
+        filename=f"{gwfname}.wel",
     )
 
     # Instantiating MODFLOW 6 output control package (flow model)
-    head_filerecord = "{}.hds".format(gwfname)
-    budget_filerecord = "{}.cbc".format(gwfname)
+    head_filerecord = f"{gwfname}.hds"
+    budget_filerecord = f"{gwfname}.cbc"
     flopy.mf6.ModflowGwfoc(
         gwf,
         head_filerecord=head_filerecord,
@@ -645,7 +645,7 @@ def build_model(sim_name, verts, cell2d, top, botm):
         sim,
         model_type="gwe6",
         modelname=gwename,
-        model_nam_file="{}.nam".format(gwename),
+        model_nam_file=f"{gwename}.nam",
     )
 
     # Create iterative model solution and register the gwe model with it
@@ -674,8 +674,8 @@ def build_model(sim_name, verts, cell2d, top, botm):
         number_orthogonalizations=2,
         preconditioner_levels=8,
         preconditioner_drop_tolerance=0.001,
-        rcloserecord="{} strict".format(rclose),
-        filename="{}.ims".format(gwename),
+        rcloserecord=f"{rclose} strict",
+        filename=f"{gwename}.ims",
     )
     sim.register_ims_package(imsgwe, [gwe.name])
 
@@ -692,18 +692,14 @@ def build_model(sim_name, verts, cell2d, top, botm):
         vertices=verts,
         cell2d=cell2d,
         pname="DISV-GWE",
-        filename="{}.disv".format(gwename),
+        filename=f"{gwename}.disv",
     )
 
     # Instantiating MODFLOW 6 heat transport initial temperature
-    flopy.mf6.ModflowGweic(
-        gwe, strt=strt_temp, pname="IC", filename="{}.ic".format(gwename)
-    )
+    flopy.mf6.ModflowGweic(gwe, strt=strt_temp, pname="IC", filename=f"{gwename}.ic")
 
     # Instantiating MODFLOW 6 heat transport advection package
-    flopy.mf6.ModflowGweadv(
-        gwe, scheme=scheme, pname="ADV", filename="{}.adv".format(gwename)
-    )
+    flopy.mf6.ModflowGweadv(gwe, scheme=scheme, pname="ADV", filename=f"{gwename}.adv")
 
     # Instantiating MODFLOW 6 heat transport energy storage package (consider renaming to est)
     prsity = determine_param(low_k_id, high_k_id, "prsity")
@@ -717,7 +713,7 @@ def build_model(sim_name, verts, cell2d, top, botm):
         heat_capacity_solid=cps,
         density_solid=rhos,
         pname="EST",
-        filename="{}.est".format(gwename),
+        filename=f"{gwename}.est",
     )
 
     # Instantiating MODFLOW 6 heat transport dispersion package
@@ -730,7 +726,7 @@ def build_model(sim_name, verts, cell2d, top, botm):
         ktw=ktw,
         kts=kts,
         pname="CND",
-        filename="{}.cnd".format(gwename),
+        filename=f"{gwename}.cnd",
     )
 
     # Instantiating MODFLOW 6 source/sink mixing package for dealing with
@@ -739,14 +735,14 @@ def build_model(sim_name, verts, cell2d, top, botm):
         ("WEL", "AUX", "TEMPERATURE"),
     ]
     flopy.mf6.ModflowGwessm(
-        gwe, sources=sourcerecarray, pname="SSM", filename="{}.ssm".format(gwename)
+        gwe, sources=sourcerecarray, pname="SSM", filename=f"{gwename}.ssm"
     )
 
     # Instantiating MODFLOW 6 heat transport output control package
     flopy.mf6.ModflowGweoc(
         gwe,
-        budget_filerecord="{}.cbc".format(gwename),
-        temperature_filerecord="{}.ucn".format(gwename),
+        budget_filerecord=f"{gwename}.cbc",
+        temperature_filerecord=f"{gwename}.ucn",
         temperatureprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
         saverecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")],
         printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")],
@@ -759,7 +755,7 @@ def build_model(sim_name, verts, cell2d, top, botm):
         exgtype="GWF6-GWE6",
         exgmnamea=gwfname,
         exgmnameb=gwename,
-        filename="{}.gwfgwe".format(gwename),
+        filename=f"{gwename}.gwfgwe",
     )
 
     return sim
@@ -950,7 +946,7 @@ def update(i):
     ani = FuncAnimation(fig, update, range(1, len(times)), init_func=init)
     # interval=25,
     writer = PillowWriter(fps=10)
-    fpth = figs_path / "{}{}".format(sim_name, ".gif")
+    fpth = figs_path / f"{sim_name}.gif"
     ani.save(fpth, writer=writer)
 
 
diff --git a/scripts/ex-gwe-barends.py b/scripts/ex-gwe-barends.py
index f15c03f3..2b4167f0 100644
--- a/scripts/ex-gwe-barends.py
+++ b/scripts/ex-gwe-barends.py
@@ -4,7 +4,7 @@
 #  Barends, F.B.J., 2010. Complete Solution for Transient Heat Transport in
 #    Porous Media, Following Lauwerier's Concept. Society of Petroleum '
 #    Engineers, Annual Technical Conference and Exhibition, Florence, Italy,
-#    19–22 September 2010.
+#    19-22 September 2010.
 #    https://doi.org/10.2118/134670-MS
 #
 # Below is a diagram of the model cell numbering (assuming delr = 4.0; values
@@ -275,7 +275,7 @@ def build_mf6_flow_model():
         top=top_overburden,
         botm=botm,
         pname="DIS",
-        filename="{}.dis".format(gwf_name),
+        filename=f"{gwf_name}.dis",
     )
 
     # Instantiate node-property flow (NPF) package
@@ -289,12 +289,12 @@ def build_mf6_flow_model():
         k=k11,
         k33=k33,
         pname="NPF-1",
-        filename="{}.npf".format(gwf_name),
+        filename=f"{gwf_name}.npf",
     )
 
     # Instantiate initial conditions package for the GWF model
     flopy.mf6.ModflowGwfic(
-        gwf, strt=top_overburden, pname="IC-1", filename="{}.ic".format(gwf_name)
+        gwf, strt=top_overburden, pname="IC-1", filename=f"{gwf_name}.ic"
     )
 
     # Instantiating MODFLOW 6 storage package
@@ -306,7 +306,7 @@ def build_mf6_flow_model():
         sy=0,
         steady_state={0: True},
         pname="STO",
-        filename="{}.sto".format(gwf_name),
+        filename=f"{gwf_name}.sto",
     )
 
     # Instantiate WEL package for adding water on the left and removing
@@ -316,19 +316,19 @@ def build_mf6_flow_model():
         auxiliary="TEMPERATURE",
         stress_period_data=wel_spd_left,
         pname="WEL-left",
-        filename="{}.wel-left".format(gwf_name),
+        filename=f"{gwf_name}.wel-left",
     )
     # Instantiate WEL package for extracting water on the right
     flopy.mf6.ModflowGwfwel(
         gwf,
         stress_period_data=wel_spd_right,
         pname="WEL-right",
-        filename="{}.wel-right".format(gwf_name),
+        filename=f"{gwf_name}.wel-right",
     )
 
     # Instantiating MODFLOW 6 output control package (flow model)
-    head_filerecord = "{}.hds".format(sim_name)
-    budget_filerecord = "{}.cbc".format(sim_name)
+    head_filerecord = f"{sim_name}.hds"
+    budget_filerecord = f"{sim_name}.cbc"
     flopy.mf6.ModflowGwfoc(
         gwf,
         head_filerecord=head_filerecord,
@@ -344,7 +344,7 @@ def build_mf6_flow_model():
 
 # +
 def build_mf6_heat_model():
-    print("Building mf6gwe model...{}".format(sim_name))
+    print(f"Building mf6gwe model...{sim_name}")
     gwename = sim_name
     sim_ws = os.path.join(workspace, sim_name, "mf6gwe")
 
@@ -355,7 +355,7 @@ def build_mf6_heat_model():
         sim,
         model_type="gwe6",
         modelname=gwename,
-        model_nam_file="{}.nam".format(gwename),
+        model_nam_file=f"{gwename}.nam",
     )
 
     # Instantiate Iterative model solution package
@@ -390,18 +390,14 @@ def build_mf6_heat_model():
         top=top_overburden,
         botm=botm,
         pname="DIS",
-        filename="{}.dis".format(gwename),
+        filename=f"{gwename}.dis",
     )
 
     # Instantiating MODFLOW 6 heat transport initial temperature
-    flopy.mf6.ModflowGweic(
-        gwe, strt=T0, pname="IC-gwe", filename="{}.ic".format(gwename)
-    )
+    flopy.mf6.ModflowGweic(gwe, strt=T0, pname="IC-gwe", filename=f"{gwename}.ic")
 
     # Instantiating MODFLOW 6 heat transport advection package
-    flopy.mf6.ModflowGweadv(
-        gwe, scheme=scheme, pname="ADV", filename="{}.adv".format(gwename)
-    )
+    flopy.mf6.ModflowGweadv(gwe, scheme=scheme, pname="ADV", filename=f"{gwename}.adv")
 
     # Instantiating MODFLOW 6 heat transport dispersion package
     if ktw != 0:
@@ -413,7 +409,7 @@ def build_mf6_heat_model():
             ktw=ktw,
             kts=kts,
             pname="CND",
-            filename="{}.cnd".format(gwename),
+            filename=f"{gwename}.cnd",
         )
 
     # Instantiating MODFLOW 6 heat transport mass storage package (consider renaming to est)
@@ -425,7 +421,7 @@ def build_mf6_heat_model():
         density_solid=rhos,
         density_water=rhow,
         pname="EST",
-        filename="{}.est".format(gwename),
+        filename=f"{gwename}.est",
     )
 
     flopy.mf6.ModflowGwectp(
@@ -435,7 +431,7 @@ def build_mf6_heat_model():
         maxbound=1,
         stress_period_data=ctp_spd,
         pname="CTP",
-        filename="{}.ctp".format(gwename),
+        filename=f"{gwename}.ctp",
     )
 
     # Instantiating MODFLOW 6 source/sink mixing package for dealing with
@@ -446,15 +442,15 @@ def build_mf6_heat_model():
         print_flows=True,
         sources=sourcerecarray,
         pname="SSM",
-        filename="{}.ssm".format(gwename),
+        filename=f"{gwename}.ssm",
     )
 
     # Instantiating MODFLOW 6 heat transport output control package
     # e.g., day 100 = 100 * 86400 = 8,640,000 = 8.64e6
     flopy.mf6.ModflowGweoc(
         gwe,
-        budget_filerecord="{}.cbc".format(gwename),
-        temperature_filerecord="{}.ucn".format(gwename),
+        budget_filerecord=f"{gwename}.cbc",
+        temperature_filerecord=f"{gwename}.ucn",
         temperatureprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
         saverecord={
             0: [
@@ -520,7 +516,7 @@ def plot_thermal_bleeding(sim_gwe):
     if plot_show:
         plt.show()
     if plot_save:
-        fpth = os.path.join(figs_path / "{}-gridView{}".format(sim_name, ".png"))
+        fpth = os.path.join(figs_path / f"{sim_name}-gridView.png")
         fig.savefig(fpth, dpi=600)
 
     # Next, plot model output
@@ -594,7 +590,7 @@ def plot_thermal_bleeding(sim_gwe):
     if plot_show:
         plt.show()
     if plot_save:
-        fpth = os.path.join(figs_path / "{}-200yrs{}".format(sim_name, ".png"))
+        fpth = os.path.join(figs_path / f"{sim_name}-200yrs.png")
         fig2.savefig(fpth, dpi=600)
 
     return
diff --git a/scripts/ex-gwe-danckwerts.py b/scripts/ex-gwe-danckwerts.py
index 1d75a51b..cd87ec82 100644
--- a/scripts/ex-gwe-danckwerts.py
+++ b/scripts/ex-gwe-danckwerts.py
@@ -366,7 +366,7 @@ def build_model(sim_name):
         scaling_method="NONE",
         reordering_method="NONE",
         relaxation_factor=relax,
-        filename="{}.ims".format(gwename),
+        filename=f"{gwename}.ims",
     )
     sim.register_ims_package(imsgwe, [gwe.name])
 
@@ -395,9 +395,7 @@ def build_model(sim_name):
     )
 
     # Instantiating MODFLOW 6 transport advection package
-    flopy.mf6.ModflowGweadv(
-        gwe, scheme=scheme, pname="ADV", filename="{}.adv".format(gwename)
-    )
+    flopy.mf6.ModflowGweadv(gwe, scheme=scheme, pname="ADV", filename=f"{gwename}.adv")
 
     # Instantiating MODFLOW 6 transport dispersion package
     flopy.mf6.ModflowGwecnd(
@@ -463,8 +461,8 @@ def build_model(sim_name):
     flopy.mf6.ModflowGweoc(
         gwe,
         pname="OC",
-        budget_filerecord="{}.cbc".format(gwename),
-        temperature_filerecord="{}.ucn".format(gwename),
+        budget_filerecord=f"{gwename}.cbc",
+        temperature_filerecord=f"{gwename}.ucn",
         temperatureprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
         saverecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")],
         printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")],
@@ -477,7 +475,7 @@ def build_model(sim_name):
         exgtype="GWF6-GWE6",
         exgmnamea=gwfname,
         exgmnameb=gwename,
-        filename="{}.gwfgwe".format(gwename),
+        filename=f"{gwename}.gwfgwe",
     )
 
     return sim
@@ -630,7 +628,7 @@ def plot_sim_vs_analytical_sln(sim):
     if plot_show:
         plt.show()
     if plot_save:
-        fpth = figs_path / "{}{}".format("ex-" + gwename + "-01", ".png")
+        fpth = figs_path / f"ex-{gwename}-01.png"
         fig.savefig(fpth, dpi=600)
 
     # Generate plots corresponding to three different times
@@ -716,7 +714,7 @@ def plot_sim_vs_analytical_sln(sim):
     if plot_show:
         plt.show()
     if plot_save:
-        fpth = figs_path / "{}{}".format("ex-" + gwename + "-02", ".png")
+        fpth = figs_path / f"ex-{gwename}-02.png"
         fig.savefig(fpth, dpi=600)
 
 
diff --git a/scripts/ex-gwe-geotherm.py b/scripts/ex-gwe-geotherm.py
index 3d0a47ba..d1694055 100644
--- a/scripts/ex-gwe-geotherm.py
+++ b/scripts/ex-gwe-geotherm.py
@@ -442,7 +442,7 @@ def build_mf6_flow_model(sim_name, silent=True):
         sim,
         modelname=gwfname,
         save_flows=True,
-        model_nam_file="{}.nam".format(gwfname),
+        model_nam_file=f"{gwfname}.nam",
     )
 
     # Instantiating MODFLOW 6 solver for flow model
@@ -471,8 +471,8 @@ def build_mf6_flow_model(sim_name, silent=True):
         number_orthogonalizations=2,
         preconditioner_levels=8,
         preconditioner_drop_tolerance=0.001,
-        rcloserecord="{} strict".format(rclose),
-        filename="{}.ims".format(sim_name),
+        rcloserecord=f"{rclose} strict",
+        filename=f"{sim_name}.ims",
     )
     sim.register_ims_package(imsgwf, [gwf.name])
 
@@ -503,7 +503,7 @@ def build_mf6_flow_model(sim_name, silent=True):
         vertices=verts,
         cell2d=cell2d,
         pname="DISV",
-        filename="{}.disv".format(gwfname),
+        filename=f"{gwfname}.disv",
     )
 
     # Instantiating MODFLOW 6 node property flow package
@@ -515,16 +515,14 @@ def build_mf6_flow_model(sim_name, silent=True):
         save_specific_discharge=True,
         save_saturation=True,
         pname="NPF",
-        filename="{}.npf".format(gwfname),
+        filename=f"{gwfname}.npf",
     )
 
     # Instantiating MODFLOW 6 initial conditions package
     flopy.mf6.ModflowGwfic(gwf, strt=strt)
 
     # Instantiating MODFLOW 6 storage package
-    flopy.mf6.ModflowGwfsto(
-        gwf, ss=0, sy=0, pname="STO", filename="{}.sto".format(gwfname)
-    )
+    flopy.mf6.ModflowGwfsto(gwf, ss=0, sy=0, pname="STO", filename=f"{gwfname}.sto")
 
     # Instantiating 1st instance of MODFLOW 6  constant head package (left side)
     # (setting auxiliary temperature to 0.0)
@@ -536,7 +534,7 @@ def build_mf6_flow_model(sim_name, silent=True):
         auxiliary="TEMPERATURE",
         stress_period_data=chd_spd,
         pname="CHD-LEFT",
-        filename="{}.left.chd".format(sim_name),
+        filename=f"{sim_name}.left.chd",
     )
 
     # Instantiating 2nd instance of MODFLOW 6 constant head package (right side)
@@ -548,12 +546,12 @@ def build_mf6_flow_model(sim_name, silent=True):
         auxiliary="TEMPERATURE",
         stress_period_data=chd_spd,
         pname="CHD-RIGHT",
-        filename="{}.right.chd".format(sim_name),
+        filename=f"{sim_name}.right.chd",
     )
 
     # Instantiating MODFLOW 6 output control package (flow model)
-    head_filerecord = "{}.hds".format(sim_name)
-    budget_filerecord = "{}.cbc".format(sim_name)
+    head_filerecord = f"{sim_name}.hds"
+    budget_filerecord = f"{sim_name}.cbc"
     flopy.mf6.ModflowGwfoc(
         gwf,
         head_filerecord=head_filerecord,
@@ -566,7 +564,7 @@ def build_mf6_flow_model(sim_name, silent=True):
 
 
 def build_mf6_heat_model(sim_name, dirichlet=0.0, neumann=0.0, silent=False):
-    print("Building mf6gwt model...{}".format(sim_name))
+    print(f"Building mf6gwt model...{sim_name}")
     gwename = "gwe-" + sim_name.split("-")[2]
     sim_ws = os.path.join(workspace, sim_name, "mf6gwe")
 
@@ -577,7 +575,7 @@ def build_mf6_heat_model(sim_name, dirichlet=0.0, neumann=0.0, silent=False):
         sim,
         model_type="gwe6",
         modelname=gwename,
-        model_nam_file="{}.nam".format(gwename),
+        model_nam_file=f"{gwename}.nam",
     )
 
     # Create iterative model solution and register the gwe model with it
@@ -606,8 +604,8 @@ def build_mf6_heat_model(sim_name, dirichlet=0.0, neumann=0.0, silent=False):
         number_orthogonalizations=2,
         preconditioner_levels=8,
         preconditioner_drop_tolerance=0.001,
-        rcloserecord="{} strict".format(rclose),
-        filename="{}.ims".format(gwename),
+        rcloserecord=f"{rclose} strict",
+        filename=f"{gwename}.ims",
     )
     sim.register_ims_package(imsgwe, [gwe.name])
 
@@ -637,14 +635,14 @@ def build_mf6_heat_model(sim_name, dirichlet=0.0, neumann=0.0, silent=False):
         vertices=verts,
         cell2d=cell2d,
         pname="DISV-GWE",
-        filename="{}.disv".format(gwename),
+        filename=f"{gwename}.disv",
     )
 
     # Instantiating MODFLOW 6 heat transport initial temperature
-    flopy.mf6.ModflowGweic(gwe, strt=strt_temp, filename="{}.ic".format(gwename))
+    flopy.mf6.ModflowGweic(gwe, strt=strt_temp, filename=f"{gwename}.ic")
 
     # Instantiating MODFLOW 6 heat transport advection package
-    flopy.mf6.ModflowGweadv(gwe, scheme=scheme, filename="{}.adv".format(gwename))
+    flopy.mf6.ModflowGweadv(gwe, scheme=scheme, filename=f"{gwename}.adv")
 
     # Instantiating MODFLOW 6 heat transport dispersion package
     if ktw != 0:
@@ -655,7 +653,7 @@ def build_mf6_heat_model(sim_name, dirichlet=0.0, neumann=0.0, silent=False):
             ktw=ktw,
             kts=kts,
             pname="CND",
-            filename="{}.dsp".format(gwename),
+            filename=f"{gwename}.dsp",
         )
 
     # Instantiating MODFLOW 6 heat transport mass storage package (consider renaming to est)
@@ -668,7 +666,7 @@ def build_mf6_heat_model(sim_name, dirichlet=0.0, neumann=0.0, silent=False):
         heat_capacity_solid=cps,
         density_solid=rhos,
         pname="EST",
-        filename="{}.est".format(gwename),
+        filename=f"{gwename}.est",
     )
 
     # Instantiating MODFLOW 6 heat transport constant temperature package (Dirichlet case)
@@ -688,7 +686,7 @@ def build_mf6_heat_model(sim_name, dirichlet=0.0, neumann=0.0, silent=False):
             stress_period_data=ctpspd,
             save_flows=False,
             pname="CTP-1",
-            filename="{}.ctp".format(gwename),
+            filename=f"{gwename}.ctp",
         )
 
     if neumann > 0:
@@ -707,7 +705,7 @@ def build_mf6_heat_model(sim_name, dirichlet=0.0, neumann=0.0, silent=False):
             stress_period_data=eslspd,
             save_flows=False,
             pname="ESR-1",
-            filename="{}.src".format(gwename),
+            filename=f"{gwename}.src",
         )
 
     # Instantiating MODFLOW 6 source/sink mixing package for dealing with
@@ -716,15 +714,13 @@ def build_mf6_heat_model(sim_name, dirichlet=0.0, neumann=0.0, silent=False):
         ("CHD-LEFT", "AUX", "TEMPERATURE"),
         ("CHD-RIGHT", "AUX", "TEMPERATURE"),
     ]
-    flopy.mf6.ModflowGwessm(
-        gwe, sources=sourcerecarray, filename="{}.ssm".format(gwename)
-    )
+    flopy.mf6.ModflowGwessm(gwe, sources=sourcerecarray, filename=f"{gwename}.ssm")
 
     # Instantiating MODFLOW 6 heat transport output control package
     flopy.mf6.ModflowGweoc(
         gwe,
-        budget_filerecord="{}.cbc".format(gwename),
-        temperature_filerecord="{}.ucn".format(gwename),
+        budget_filerecord=f"{gwename}.cbc",
+        temperature_filerecord=f"{gwename}.ucn",
         temperatureprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
         saverecord={
             49: [("TEMPERATURE", "LAST"), ("BUDGET", "LAST")],
@@ -797,7 +793,7 @@ def plot_grid(sim):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = os.path.join(figs_path / "{}-grid{}".format(simname, ".png"))
+            fpth = os.path.join(figs_path / f"{simname}-grid.png")
             fig.savefig(fpth, dpi=300)
 
 
@@ -828,7 +824,7 @@ def plot_grid_inset(sim):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = os.path.join(figs_path / "{}-grid-inset{}".format(simname, ".png"))
+            fpth = os.path.join(figs_path / f"{simname}-grid-inset.png")
             fig.savefig(fpth, dpi=300)
 
 
@@ -857,7 +853,7 @@ def plot_head(sim):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = os.path.join(figs_path / "{}-head{}".format(simname, ".png"))
+            fpth = os.path.join(figs_path / f"{simname}-head.png")
             fig.savefig(fpth, dpi=300)
 
 
@@ -951,7 +947,7 @@ def plot_temperature(sim, scen, time_):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = os.path.join(figs_path / "{}-temp50days{}".format(simname, ".png"))
+            fpth = os.path.join(figs_path / f"{simname}-temp50days.png")
             fig.savefig(fpth, dpi=300)
 
 
diff --git a/scripts/ex-gwe-prt.py b/scripts/ex-gwe-prt.py
index dcc2a887..895d876a 100644
--- a/scripts/ex-gwe-prt.py
+++ b/scripts/ex-gwe-prt.py
@@ -208,10 +208,10 @@ def build_gwf_sim(name):
         sy=0,
         steady_state={0: True},
         pname="STO",
-        filename="{}.sto".format(gwf_name),
+        filename=f"{gwf_name}.sto",
     )
 
-    flopy.mf6.ModflowGwfic(gwf, strt=1.0, pname="IC", filename="{}.ic".format(gwf_name))
+    flopy.mf6.ModflowGwfic(gwf, strt=1.0, pname="IC", filename=f"{gwf_name}.ic")
 
     flopy.mf6.ModflowGwfoc(
         gwf,
@@ -219,7 +219,7 @@ def build_gwf_sim(name):
         head_filerecord=f"{gwf_name}.hds",
         saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
         printrecord=[("HEAD", "LAST"), ("BUDGET", "LAST")],
-        filename="{}.oc".format(gwf_name),
+        filename=f"{gwf_name}.oc",
     )
 
     chdlist = []
@@ -241,7 +241,7 @@ def build_gwf_sim(name):
         auxiliary="TEMPERATURE",
         stress_period_data=chdlist,
         pname="CHD",
-        filename="{}.chd".format(gwf_name),
+        filename=f"{gwf_name}.chd",
     )
 
     return sim
@@ -258,7 +258,7 @@ def build_gwe_sim(name):
         sim,
         model_type="gwe6",
         modelname=gwe_name,
-        model_nam_file="{}.nam".format(gwe_name),
+        model_nam_file=f"{gwe_name}.nam",
     )
     imsgwe = flopy.mf6.ModflowIms(
         sim,
@@ -284,14 +284,10 @@ def build_gwe_sim(name):
         top=top,
         botm=botm,
         pname="DISV",
-        filename="{}.disv".format(gwe_name),
-    )
-    flopy.mf6.ModflowGweic(
-        gwe, strt=strt_temp, pname="IC", filename="{}.ic".format(gwe_name)
-    )
-    flopy.mf6.ModflowGweadv(
-        gwe, scheme=scheme, pname="ADV", filename="{}.adv".format(gwe_name)
+        filename=f"{gwe_name}.disv",
     )
+    flopy.mf6.ModflowGweic(gwe, strt=strt_temp, pname="IC", filename=f"{gwe_name}.ic")
+    flopy.mf6.ModflowGweadv(gwe, scheme=scheme, pname="ADV", filename=f"{gwe_name}.adv")
     if ktw != 0:
         flopy.mf6.ModflowGwecnd(
             gwe,
@@ -300,7 +296,7 @@ def build_gwe_sim(name):
             ktw=ktw,
             kts=kts,
             pname="CND-e",
-            filename="{}.cnd".format(gwe_name),
+            filename=f"{gwe_name}.cnd",
         )
 
     flopy.mf6.ModflowGweest(
@@ -312,7 +308,7 @@ def build_gwe_sim(name):
         heat_capacity_solid=cps,
         density_solid=rhos,
         pname="EST-e",
-        filename="{}.est".format(gwe_name),
+        filename=f"{gwe_name}.est",
     )
 
     sourcerecarray = [
@@ -320,12 +316,12 @@ def build_gwe_sim(name):
         ("CHD", "AUX", "TEMPERATURE"),
     ]
     flopy.mf6.ModflowGwessm(
-        gwe, sources=sourcerecarray, pname="SSM-e", filename="{}.ssm".format(gwe_name)
+        gwe, sources=sourcerecarray, pname="SSM-e", filename=f"{gwe_name}.ssm"
     )
     flopy.mf6.ModflowGweoc(
         gwe,
-        budget_filerecord="{}.cbc".format(gwe_name),
-        temperature_filerecord="{}.ucn".format(gwe_name),
+        budget_filerecord=f"{gwe_name}.cbc",
+        temperature_filerecord=f"{gwe_name}.ucn",
         temperatureprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
         saverecord={0: [("TEMPERATURE", "ALL"), ("BUDGET", "ALL")]},
         printrecord=[("TEMPERATURE", "LAST"), ("BUDGET", "LAST")],
@@ -635,7 +631,7 @@ def plot_results(gwf_sim, gwe_sim, prt_sim, silent=True):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(sim_name, ".png")
+            fpth = figs_path / f"{sim_name}.png"
             fig.savefig(fpth, dpi=600)
 
 
diff --git a/scripts/ex-gwe-radial.py b/scripts/ex-gwe-radial.py
index 512c4757..97c52dbe 100644
--- a/scripts/ex-gwe-radial.py
+++ b/scripts/ex-gwe-radial.py
@@ -537,7 +537,7 @@ def build_mf6_flow_model(sim_name, left_chd_spd=None, right_chd_spd=None, silent
         sim,
         modelname=gwfname,
         save_flows=True,
-        model_nam_file="{}.nam".format(gwfname),
+        model_nam_file=f"{gwfname}.nam",
     )
 
     # Instantiating MODFLOW 6 solver for flow model
@@ -566,8 +566,8 @@ def build_mf6_flow_model(sim_name, left_chd_spd=None, right_chd_spd=None, silent
         number_orthogonalizations=2,
         preconditioner_levels=8,
         preconditioner_drop_tolerance=0.001,
-        rcloserecord="{} strict".format(rclose),
-        filename="{}.ims".format(sim_name),
+        rcloserecord=f"{rclose} strict",
+        filename=f"{sim_name}.ims",
     )
     sim.register_ims_package(imsgwf, [gwf.name])
 
@@ -585,7 +585,7 @@ def build_mf6_flow_model(sim_name, left_chd_spd=None, right_chd_spd=None, silent
         vertices=verts,
         cell2d=cell2d,
         pname="DISV",
-        filename="{}.disv".format(gwfname),
+        filename=f"{gwfname}.disv",
     )
 
     # Instantiating MODFLOW 6 node property flow package
@@ -597,7 +597,7 @@ def build_mf6_flow_model(sim_name, left_chd_spd=None, right_chd_spd=None, silent
         save_specific_discharge=True,
         save_saturation=True,
         pname="NPF",
-        filename="{}.npf".format(gwfname),
+        filename=f"{gwfname}.npf",
     )
 
     # Instantiating MODFLOW 6 initial conditions package
@@ -612,7 +612,7 @@ def build_mf6_flow_model(sim_name, left_chd_spd=None, right_chd_spd=None, silent
         sy=0,
         steady_state={0: True},
         pname="STO",
-        filename="{}.sto".format(gwfname),
+        filename=f"{gwfname}.sto",
     )
 
     # Instantiating 1st instance of MODFLOW 6 constant head package (left side)
@@ -623,7 +623,7 @@ def build_mf6_flow_model(sim_name, left_chd_spd=None, right_chd_spd=None, silent
         auxiliary="TEMPERATURE",
         stress_period_data=left_chd_spd,
         pname="CHD-LEFT",
-        filename="{}.left.chd".format(sim_name),
+        filename=f"{sim_name}.left.chd",
     )
 
     right_chd_spd = {0: right_chd_spd}
@@ -632,12 +632,12 @@ def build_mf6_flow_model(sim_name, left_chd_spd=None, right_chd_spd=None, silent
         auxiliary="TEMPERATURE",
         stress_period_data=right_chd_spd,
         pname="CHD-RIGHT",
-        filename="{}.right.chd".format(sim_name),
+        filename=f"{sim_name}.right.chd",
     )
 
     # Instantiating MODFLOW 6 output control package (flow model)
-    head_filerecord = "{}.hds".format(sim_name)
-    budget_filerecord = "{}.cbc".format(sim_name)
+    head_filerecord = f"{sim_name}.hds"
+    budget_filerecord = f"{sim_name}.cbc"
     flopy.mf6.ModflowGwfoc(
         gwf,
         head_filerecord=head_filerecord,
@@ -656,7 +656,7 @@ def build_mf6_heat_model(
     gwvelocity=0.0,
     silent=False,
 ):
-    print("Building mf6gwt model...{}".format(sim_name))
+    print(f"Building mf6gwt model...{sim_name}")
     gwename = "gwe-" + sim_name.split("-")[2]
     sim_ws = os.path.join(workspace, sim_name[:-2], "mf6gwe" + scen_ext)
     sim = flopy.mf6.MFSimulation(sim_name=sim_name, sim_ws=sim_ws, exe_name="mf6")
@@ -666,7 +666,7 @@ def build_mf6_heat_model(
         sim,
         model_type="gwe6",
         modelname=gwename,
-        model_nam_file="{}.nam".format(gwename),
+        model_nam_file=f"{gwename}.nam",
     )
 
     # Create iterative model solution and register the gwe model with it
@@ -695,8 +695,8 @@ def build_mf6_heat_model(
         number_orthogonalizations=2,
         preconditioner_levels=8,
         preconditioner_drop_tolerance=0.001,
-        rcloserecord="{} strict".format(rclose),
-        filename="{}.ims".format(gwename),
+        rcloserecord=f"{rclose} strict",
+        filename=f"{gwename}.ims",
     )
     sim.register_ims_package(imsgwe, [gwe.name])
 
@@ -731,14 +731,14 @@ def build_mf6_heat_model(
         vertices=verts,
         cell2d=cell2d,
         pname="DISV-GWE",
-        filename="{}.disv".format(gwename),
+        filename=f"{gwename}.disv",
     )
 
     # Instantiating MODFLOW 6 heat transport initial temperature
-    flopy.mf6.ModflowGweic(gwe, strt=strt_temp, filename="{}.ic".format(gwename))
+    flopy.mf6.ModflowGweic(gwe, strt=strt_temp, filename=f"{gwename}.ic")
 
     # Instantiating MODFLOW 6 heat transport advection package
-    flopy.mf6.ModflowGweadv(gwe, scheme=scheme, filename="{}.adv".format(gwename))
+    flopy.mf6.ModflowGweadv(gwe, scheme=scheme, filename=f"{gwename}.adv")
 
     # Instantiating MODFLOW 6 heat transport dispersion package
     if ktw != 0:
@@ -749,7 +749,7 @@ def build_mf6_heat_model(
             ktw=ktw,
             kts=kts,
             pname="CND",
-            filename="{}.cnd".format(gwename),
+            filename=f"{gwename}.cnd",
         )
 
     # Instantiating MODFLOW 6 heat transport mass storage package (consider renaming to est)
@@ -762,7 +762,7 @@ def build_mf6_heat_model(
         heat_capacity_solid=cps,
         density_solid=rhos,
         pname="EST",
-        filename="{}.est".format(gwename),
+        filename=f"{gwename}.est",
     )
 
     # Instantiating MODFLOW 6 heat transport constant temperature package (Dirichlet case)
@@ -777,7 +777,7 @@ def build_mf6_heat_model(
             stress_period_data=ctpspd,
             save_flows=False,
             pname="CTP-1",
-            filename="{}.ctp".format(gwename),
+            filename=f"{gwename}.ctp",
         )
 
     if neumann > 0 and dirichlet == 0:
@@ -791,7 +791,7 @@ def build_mf6_heat_model(
             stress_period_data=esl_spd,
             save_flows=False,
             pname="ESL-1",
-            filename="{}.esl".format(gwename),
+            filename=f"{gwename}.esl",
         )
 
     # Instantiating MODFLOW 6 source/sink mixing package for dealing with
@@ -800,15 +800,13 @@ def build_mf6_heat_model(
         ("CHD-LEFT", "AUX", "TEMPERATURE"),
         ("CHD-RIGHT", "AUX", "TEMPERATURE"),
     ]
-    flopy.mf6.ModflowGwessm(
-        gwe, sources=sourcerecarray, filename="{}.ssm".format(gwename)
-    )
+    flopy.mf6.ModflowGwessm(gwe, sources=sourcerecarray, filename=f"{gwename}.ssm")
 
     # Instantiating MODFLOW 6 heat transport output control package
     flopy.mf6.ModflowGweoc(
         gwe,
-        budget_filerecord="{}.cbc".format(gwename),
-        temperature_filerecord="{}.ucn".format(gwename),
+        budget_filerecord=f"{gwename}.cbc",
+        temperature_filerecord=f"{gwename}.ucn",
         temperatureprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
         saverecord={47: [("TEMPERATURE", "LAST"), ("BUDGET", "LAST")]},
         printrecord=[("TEMPERATURE", "LAST"), ("BUDGET", "LAST")],
@@ -875,7 +873,7 @@ def plot_grid(sim):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}-grid{}".format(simname, ".png")
+            fpth = figs_path / f"{simname}-grid.png"
             fig.savefig(fpth, dpi=300)
 
 
@@ -907,7 +905,7 @@ def plot_grid_inset(sim):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}-gridinset{}".format(simname, ".png")
+            fpth = figs_path / f"{simname}-gridinset.png"
             fig.savefig(fpth, dpi=300)
 
 
@@ -935,7 +933,7 @@ def plot_head(sim):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}-head{}".format(simname, ".png")
+            fpth = figs_path / f"{simname}-head.png"
             fig.savefig(fpth, dpi=300)
 
 
@@ -1010,7 +1008,7 @@ def plot_temperature(sim, idx, scen_txt, vel_txt):
             plt.show()
         if plot_save:
             fname = sim.name
-            fpth = figs_path / "{}-temp48{}".format(sim_name, ".png")
+            fpth = figs_path / f"{sim_name}-temp48.png"
             fig.savefig(fpth, dpi=300)
 
 
diff --git a/scripts/ex-gwe-vsc.py b/scripts/ex-gwe-vsc.py
index 14c05cba..b30424be 100644
--- a/scripts/ex-gwe-vsc.py
+++ b/scripts/ex-gwe-vsc.py
@@ -155,7 +155,7 @@
 # A function for customizing colormaps
 def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100):
     new_cmap = colors.LinearSegmentedColormap.from_list(
-        "trunc({n},{a:.2f},{b:.2f})".format(n=cmap.name, a=minval, b=maxval),
+        f"trunc({cmap.name},{minval:.2f},{maxval:.2f})",
         cmap(np.linspace(minval, maxval, n)),
     )
     return new_cmap
@@ -287,7 +287,7 @@ def add_gwf_model(sim, gwfname, activate_vsc):
 
     if activate_vsc:
         # Instantiate viscosity (VSC) package
-        vsc_filerecord = "{}.vsc.bin".format(gwfname)
+        vsc_filerecord = f"{gwfname}.vsc.bin"
         gwename2 = sim_name + "-02"
         vsc_pd = [(iviscspec, dviscdc, cviscref, gwename2, "TEMPERATURE")]
         flopy.mf6.ModflowGwfvsc(
@@ -301,7 +301,7 @@ def add_gwf_model(sim, gwfname, activate_vsc):
             nviscspecies=len(vsc_pd),
             packagedata=vsc_pd,
             pname="VSC",
-            filename="{}.vsc".format(gwfname),
+            filename=f"{gwfname}.vsc",
         )
 
     # Instantiating MODFLOW 6 output control package for flow model
@@ -456,7 +456,7 @@ def add_gwe_model(sim, gwename, temp_upper=4.0, temp_lower=4.0):
         scaling_method="NONE",
         reordering_method="NONE",
         relaxation_factor=relax,
-        filename="{}.ims".format(gwename),
+        filename=f"{gwename}.ims",
     )
     sim.register_ims_package(imsgwe, [gwename])
 
@@ -473,7 +473,7 @@ def add_gwe_model(sim, gwename, temp_upper=4.0, temp_lower=4.0):
         botm=botm,
         idomain=1,
         pname="DIS",
-        filename="{}.dis".format(gwename),
+        filename=f"{gwename}.dis",
     )
 
     # Instantiate Mobile Storage and Transfer package
@@ -487,7 +487,7 @@ def add_gwe_model(sim, gwename, temp_upper=4.0, temp_lower=4.0):
         heat_capacity_solid=cps,
         density_solid=rhos,
         pname="EST",
-        filename="{}.est".format(gwename),
+        filename=f"{gwename}.est",
     )
 
     # Instantiate Energy Transport Initial Conditions package
@@ -502,7 +502,7 @@ def add_gwe_model(sim, gwename, temp_upper=4.0, temp_lower=4.0):
         maxbound=len(ctp_spd),
         stress_period_data=ctp_spd,
         pname="CTP",
-        filename="{}.ctp".format(gwename),
+        filename=f"{gwename}.ctp",
     )
 
     # Instantiate Dispersion package (also handles conduction)
@@ -512,7 +512,7 @@ def add_gwe_model(sim, gwename, temp_upper=4.0, temp_lower=4.0):
         ktw=ktw,
         kts=kts,
         pname="CND",
-        filename="{}.cnd".format(gwename),
+        filename=f"{gwename}.cnd",
     )
 
     # Instantiating MODFLOW 6 transport source-sink mixing package
@@ -522,17 +522,17 @@ def add_gwe_model(sim, gwename, temp_upper=4.0, temp_lower=4.0):
         ("WEL-1", "AUX", "TEMPERATURE"),
     ]
     flopy.mf6.ModflowGwessm(
-        gwe, sources=sourcerecarray, pname="SSM", filename="{}.ssm".format(gwename)
+        gwe, sources=sourcerecarray, pname="SSM", filename=f"{gwename}.ssm"
     )
 
     # Instantiate Output Control package for transport
     flopy.mf6.ModflowGweoc(
         gwe,
-        temperature_filerecord="{}.ucn".format(gwename),
+        temperature_filerecord=f"{gwename}.ucn",
         saverecord=[("TEMPERATURE", "ALL")],
         temperatureprintrecord=[("COLUMNS", 3, "WIDTH", 20, "DIGITS", 8, "GENERAL")],
         printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")],
-        filename="{}.oc".format(gwename),
+        filename=f"{gwename}.oc",
     )
 
     return sim
@@ -589,7 +589,7 @@ def build_models(sim_name, silent=False, temp_upper=4.0, temp_lower=4.0):
                 exgtype="GWF6-GWE6",
                 exgmnamea=gwfname,
                 exgmnameb=gwename,
-                filename="{}.gwfgwe".format(gwename),
+                filename=f"{gwename}.gwfgwe",
             )
 
     return sim
diff --git a/scripts/ex-gwf-advtidal.py b/scripts/ex-gwf-advtidal.py
index 25391e29..56120fc6 100644
--- a/scripts/ex-gwf-advtidal.py
+++ b/scripts/ex-gwf-advtidal.py
@@ -3,7 +3,7 @@
 # This model represents highlands bordered by a tidal estuary.  The model
 # has 3 layers representing an upper aquifer, confining bed, and lower aquifer.
 # The grid is 15 rows by 10 columns.  The length unit is feet and the time unit
-# is days. Each cell is 500 ft × 500 ft.  The estuary is represented by GHB
+# is days. Each cell is 500 ft x 500 ft.  The estuary is represented by GHB
 # boundaries in column 10.  Two rivers cross the area from left to right.
 # Recharge is zoned by the use of three Recharge-Package input files
 
@@ -452,7 +452,7 @@ def plot_ts(sim):
             if plot_show:
                 plt.show()
             if plot_save:
-                fpth = figs_path / "{}-{}{}".format(sim_name, obs_fig[iplot], ".png")
+                fpth = figs_path / f"{sim_name}-{obs_fig[iplot]}.png"
                 fig.savefig(fpth)
 
 
diff --git a/scripts/ex-gwf-bcf2ss.py b/scripts/ex-gwf-bcf2ss.py
index 993938a6..a30144e8 100644
--- a/scripts/ex-gwf-bcf2ss.py
+++ b/scripts/ex-gwf-bcf2ss.py
@@ -346,7 +346,7 @@ def plot_results(silent=True):
         verbosity_level = 1
 
     with styles.USGSMap():
-        name = list(parameters.keys())[0]
+        name = next(iter(parameters.keys()))
         sim_ws = os.path.join(workspace, name)
         sim = flopy.mf6.MFSimulation.load(
             sim_name=sim_name, sim_ws=sim_ws, verbosity_level=verbosity_level
diff --git a/scripts/ex-gwf-csub-p02.py b/scripts/ex-gwf-csub-p02.py
index 22aaa59e..fdb123e4 100644
--- a/scripts/ex-gwf-csub-p02.py
+++ b/scripts/ex-gwf-csub-p02.py
@@ -560,14 +560,14 @@ def plot_comp_q_comparison(sim, silent=True):
                 ax.set_xlim(0, 100)
                 ylabel = (
                     "Head-based minus effective stress-based\nsubsidence, "
-                    + "in % of ultimate value"
+                    "in % of ultimate value"
                 )
             else:
                 ax.set_ylim(0, 8)
                 ax.set_xlim(0, 100)
                 ylabel = (
                     "Top minus bottom interbed effective stress-\nbased "
-                    + "rate, in % of head-based drainage rate"
+                    "rate, in % of head-based drainage rate"
                 )
             ax.set_xlabel("Percent of time constant")
             ax.set_ylabel(ylabel)
@@ -652,7 +652,7 @@ def plot_head_comparison(sim, silent=True):
                 styles.heading(ax, letter=chr(ord("A") + idx))
                 ax.set_yticks(yticks)
                 styles.remove_edge_ticks(ax)
-                text = r"$\frac{t}{\tau_0}$ = " + "{}".format(pct_vals[idx] / 100.0)
+                text = r"$\frac{t}{\tau_0}$ = " + f"{pct_vals[idx] / 100.0}"
                 ax.text(
                     0.25,
                     0.01,
@@ -671,8 +671,8 @@ def plot_head_comparison(sim, silent=True):
             else:
                 if idx == 2:
                     text = (
-                        "Difference in head-based and effective stress-based"
-                        + "\ninterbed heads, in percent of head-based interbed heads"
+                        "Difference in head-based and effective stress-based\n"
+                        "interbed heads, in percent of head-based interbed heads"
                     )
                     ax.set_xlabel(text)
                 ax.set_yticklabels([])
diff --git a/scripts/ex-gwf-csub-p03.py b/scripts/ex-gwf-csub-p03.py
index 18672f5d..548c9559 100644
--- a/scripts/ex-gwf-csub-p03.py
+++ b/scripts/ex-gwf-csub-p03.py
@@ -679,7 +679,7 @@ def export_tables(silent=True):
         col_widths = (0.1, 0.25)
         build_table(caption, fpth, arr, headings=headings, col_widths=col_widths)
 
-        caption = "Interbed storage properties for example {}.".format(sim_name)
+        caption = f"Interbed storage properties for example {sim_name}."
         headings = (
             "Interbed",
             "Layer",
@@ -1023,7 +1023,7 @@ def process_dtw_obs(fpth):
                 v[key] *= -1.0
             return v
 
-        name = list(parameters.keys())[0]
+        name = next(iter(parameters.keys()))
         pth = os.path.join(workspace, name, f"{name}.gwf.obs.csv")
         hdata = process_dtw_obs(pth)
 
@@ -1067,7 +1067,7 @@ def process_dtw_obs(fpth):
 
 def plot_head_es_comparison(silent=True):
     with styles.USGSPlot() as fs:
-        name = list(parameters.keys())[0]
+        name = next(iter(parameters.keys()))
         pth = os.path.join(workspace, name, f"{name}.csub.obs.csv")
         hb = process_csub_obs(pth)
 
@@ -1362,7 +1362,7 @@ def plot_calibration(silent=True):
         )
 
         for idx, ixc in enumerate(ixs):
-            text = "{}".format(df_iobs_pc.index[ixc].strftime("%m/%d/%Y"))
+            text = f"{df_iobs_pc.index[ixc]:%m/%d/%Y}"
             if df_iobs_pc.index[ixc].month == 4:
                 dxc = -0.001
                 dyc = -1
@@ -1392,8 +1392,8 @@ def plot_calibration(silent=True):
             df_sim_pc.index[0].strftime("%B %d, %Y"), length_units
         )
         ytext = (
-            "Effective stress at the bottom of\nthe lower aquifer, in "
-            + f"{length_units} of water"
+            "Effective stress at the bottom of\n"
+            f"the lower aquifer, in {length_units} of water"
         )
         ax.set_xlabel(xtext)
         ax.set_ylabel(ytext)
diff --git a/scripts/ex-gwf-curvilinear-90.py b/scripts/ex-gwf-curvilinear-90.py
index 79c13dcf..b1b5d467 100644
--- a/scripts/ex-gwf-curvilinear-90.py
+++ b/scripts/ex-gwf-curvilinear-90.py
@@ -24,7 +24,6 @@
 import pathlib as pl
 from itertools import cycle
 from math import sqrt
-from typing import Dict, List
 
 import flopy
 import git
@@ -153,10 +152,10 @@ class DisvPropertyContainer:
     nlay: int
     ncpl: int
     nvert: int
-    vertices: List[list]  # [[iv, xv, yv], ...]
-    cell2d: List[list]  # [[ic, xc, yc, ncvert, icvert], ...]
+    vertices: list[list]  # [[iv, xv, yv], ...]
+    cell2d: list[list]  # [[ic, xc, yc, ncvert, icvert], ...]
     top: np.ndarray
-    botm: List[np.ndarray]
+    botm: list[np.ndarray]
     origin_x: float
     origin_y: float
     rotation: float
@@ -1166,7 +1165,7 @@ class DisvGridMerger:
         kwargs to DisvPropertyContainer.plot_grid(...).
     """
 
-    grids: Dict[str, DisvPropertyContainer]
+    grids: dict[str, DisvPropertyContainer]
     merged: DisvPropertyContainer
     snap_vertices: dict
     connect_tolerance: dict
@@ -1333,9 +1332,7 @@ def _get_vertex_xy(self, name, iv):
         for iv_orig, xv, yv in vertices:
             if iv == iv_orig:
                 return xv, yv
-        raise RuntimeError(
-            "DisvGridMerger: " f"Failed to find vertex {iv} in grid {name}"
-        )
+        raise RuntimeError(f"DisvGridMerger: Failed to find vertex {iv} in grid {name}")
 
     def _find_merged_vertex(self, xv, yv, tol):
         for iv, xv_chk, yv_chk in self.merged.vertices:
@@ -1350,7 +1347,7 @@ def _replace_vertex_xy(self, iv, xv, yv):
                 vert[2] = yv
                 return
         raise RuntimeError(
-            "DisvGridMerger: Unknown code error - " f"failed to locate vertex {iv}"
+            f"DisvGridMerger: Unknown code error - failed to locate vertex {iv}"
         )
 
     def _clear_attribute(self):
@@ -2459,7 +2456,7 @@ def plot_analytical(sim, verbose=False):
         styles.graph_legend(ax)
         fig.tight_layout()
         if plot_save:
-            fpth = figs_path / "{}-{}{}".format(sim_name, obs_fig, ".png")
+            fpth = figs_path / f"{sim_name}-{obs_fig}.png"
             fig.savefig(fpth)
 
 
diff --git a/scripts/ex-gwf-curvilinear.py b/scripts/ex-gwf-curvilinear.py
index f09a856d..6bcd9fdf 100644
--- a/scripts/ex-gwf-curvilinear.py
+++ b/scripts/ex-gwf-curvilinear.py
@@ -24,7 +24,6 @@
 import pathlib as pl
 from itertools import cycle
 from math import sqrt
-from typing import Dict, List
 
 import flopy
 import git
@@ -153,10 +152,10 @@ class DisvPropertyContainer:
     nlay: int
     ncpl: int
     nvert: int
-    vertices: List[list]  # [[iv, xv, yv], ...]
-    cell2d: List[list]  # [[ic, xc, yc, ncvert, icvert], ...]
+    vertices: list[list]  # [[iv, xv, yv], ...]
+    cell2d: list[list]  # [[ic, xc, yc, ncvert, icvert], ...]
     top: np.ndarray
-    botm: List[np.ndarray]
+    botm: list[np.ndarray]
     origin_x: float
     origin_y: float
     rotation: float
@@ -1166,7 +1165,7 @@ class DisvGridMerger:
         kwargs to DisvPropertyContainer.plot_grid(...).
     """
 
-    grids: Dict[str, DisvPropertyContainer]
+    grids: dict[str, DisvPropertyContainer]
     merged: DisvPropertyContainer
     snap_vertices: dict
     connect_tolerance: dict
@@ -1333,9 +1332,7 @@ def _get_vertex_xy(self, name, iv):
         for iv_orig, xv, yv in vertices:
             if iv == iv_orig:
                 return xv, yv
-        raise RuntimeError(
-            "DisvGridMerger: " f"Failed to find vertex {iv} in grid {name}"
-        )
+        raise RuntimeError(f"DisvGridMerger: Failed to find vertex {iv} in grid {name}")
 
     def _find_merged_vertex(self, xv, yv, tol):
         for iv, xv_chk, yv_chk in self.merged.vertices:
@@ -1350,7 +1347,7 @@ def _replace_vertex_xy(self, iv, xv, yv):
                 vert[2] = yv
                 return
         raise RuntimeError(
-            "DisvGridMerger: Unknown code error - " f"failed to locate vertex {iv}"
+            f"DisvGridMerger: Unknown code error - failed to locate vertex {iv}"
         )
 
     def _clear_attribute(self):
diff --git a/scripts/ex-gwf-drn-p01.py b/scripts/ex-gwf-drn-p01.py
index 4c8d66f1..20b6a596 100644
--- a/scripts/ex-gwf-drn-p01.py
+++ b/scripts/ex-gwf-drn-p01.py
@@ -589,7 +589,7 @@ def plot_gwseep_results(silent=True):
     """Plot groundwater seepage results"""
     with styles.USGSPlot() as fs:
         # load the observations
-        name = list(parameters.keys())[0]
+        name = next(iter(parameters.keys()))
         fpth = os.path.join(workspace, name, f"{sim_name}.surfrate.obs.csv")
         drn = flopy.utils.Mf6Obs(fpth).data
         name = list(parameters.keys())[1]
@@ -698,7 +698,7 @@ def plot_gwseep_results(silent=True):
 
 def export_tables(silent=True):
     if plot_save:
-        caption = "Infiltration and pumping rates for example {}.".format(sim_name)
+        caption = f"Infiltration and pumping rates for example {sim_name}."
         headings = (
             "Stress period",
             "Infiltration rate",
diff --git a/scripts/ex-gwf-fhb.py b/scripts/ex-gwf-fhb.py
index 9defcba2..c3cb4535 100644
--- a/scripts/ex-gwf-fhb.py
+++ b/scripts/ex-gwf-fhb.py
@@ -251,7 +251,7 @@ def plot_ts(sim):
             ax.set_ylabel(ylabel[iplot])
             styles.graph_legend(ax)
             if plot_save:
-                fpth = figs_path / "{}-{}{}".format(sim_name, obs_fig[iplot], ".png")
+                fpth = figs_path / f"{sim_name}-{obs_fig[iplot]}.png"
                 fig.savefig(fpth)
 
 
diff --git a/scripts/ex-gwf-lgr.py b/scripts/ex-gwf-lgr.py
index a744b9f4..e88d043b 100644
--- a/scripts/ex-gwf-lgr.py
+++ b/scripts/ex-gwf-lgr.py
@@ -861,7 +861,7 @@ def plot_results(mf6, idx):
     with styles.USGSPlot():
         # Start by retrieving some output
         mf6_out_pth = mf6.simulation_data.mfpath.get_sim_path()
-        sfr_parent_bud_file = list(mf6.model_names)[0] + ".sfr.bud"
+        sfr_parent_bud_file = next(iter(mf6.model_names)) + ".sfr.bud"
         sfr_child_bud_file = list(mf6.model_names)[1] + ".sfr.bud"
         sfr_parent_out = os.path.join(mf6_out_pth, sfr_parent_bud_file)
         sfr_child_out = os.path.join(mf6_out_pth, sfr_child_bud_file)
@@ -920,7 +920,7 @@ def plot_results(mf6, idx):
             strmQ[i + (j + 1)] = qc[0][0][jtm][2]
 
         # The flow that is passed between the parent and child grids comes next
-        from_mvr = [itm[2] for itm in dat_fmp[0] if itm[2] > 0][0]
+        from_mvr = next(itm[2] for itm in dat_fmp[0] if itm[2] > 0)
         strmQ[i + j + 2] = from_mvr
 
         # Finally, process the remaining parent model stream reaches
diff --git a/scripts/ex-gwf-maw-p01.py b/scripts/ex-gwf-maw-p01.py
index 42596950..a71ea7ee 100644
--- a/scripts/ex-gwf-maw-p01.py
+++ b/scripts/ex-gwf-maw-p01.py
@@ -210,7 +210,7 @@ def run_models(sim, silent=True):
 def plot_maw_results(silent=True):
     with styles.USGSPlot():
         # load the observations
-        name = list(parameters.keys())[0]
+        name = next(iter(parameters.keys()))
         fpth = os.path.join(workspace, name, f"{sim_name}.maw.obs.csv")
         maw0 = flopy.utils.Mf6Obs(fpth).data
         name = list(parameters.keys())[1]
@@ -287,7 +287,7 @@ def plot_grid(silent=True):
         verbosity = 0
     else:
         verbosity = 1
-    name = list(parameters.keys())[0]
+    name = next(iter(parameters.keys()))
     sim_ws = os.path.join(workspace, name)
     sim = flopy.mf6.MFSimulation.load(
         sim_name=sim_name, sim_ws=sim_ws, verbosity_level=verbosity
diff --git a/scripts/ex-gwf-maw-p03.py b/scripts/ex-gwf-maw-p03.py
index 4e303a7b..bff62e0a 100644
--- a/scripts/ex-gwf-maw-p03.py
+++ b/scripts/ex-gwf-maw-p03.py
@@ -281,7 +281,7 @@ def build_regional(name):
 
 def build_local(name, simulation):
     # get regional heads for constant head boundaries
-    pth = list(parameters.keys())[0]
+    pth = next(iter(parameters.keys()))
     fpth = os.path.join(workspace, pth, f"{sim_name}.hds")
     try:
         h = flopy.utils.HeadFile(fpth).get_data()
@@ -440,7 +440,7 @@ def plot_maw_results(silent=True):
         hgwf /= ihds
 
         if silent:
-            print("MAW head: {}  Average head: {}".format(maw["H0"], hgwf))
+            print(f"MAW head: {maw['H0']}  Average head: {hgwf}")
 
         zelev = sorted(list(set(list(obs_elev.values()))), reverse=True)
 
@@ -529,7 +529,7 @@ def plot_regional_grid(silent=True):
         verbosity = 0
     else:
         verbosity = 1
-    name = list(parameters.keys())[0]
+    name = next(iter(parameters.keys()))
     sim_ws = os.path.join(workspace, name)
     sim = flopy.mf6.MFSimulation.load(
         sim_name=sim_name, sim_ws=sim_ws, verbosity_level=verbosity
diff --git a/scripts/ex-gwf-nwt-p02.py b/scripts/ex-gwf-nwt-p02.py
index 7354f4bf..c4dcc830 100644
--- a/scripts/ex-gwf-nwt-p02.py
+++ b/scripts/ex-gwf-nwt-p02.py
@@ -243,7 +243,7 @@ def plot_results(silent=True):
 
     with styles.USGSMap():
         # load the newton model
-        name = list(parameters.keys())[0]
+        name = next(iter(parameters.keys()))
         sim_ws = os.path.join(workspace, name)
         sim = flopy.mf6.MFSimulation.load(
             sim_name=sim_name, sim_ws=sim_ws, verbosity_level=verbosity_level
diff --git a/scripts/ex-gwf-radial.py b/scripts/ex-gwf-radial.py
index 8f9e3b3c..c6b5c569 100644
--- a/scripts/ex-gwf-radial.py
+++ b/scripts/ex-gwf-radial.py
@@ -370,8 +370,8 @@ def __init__(
         if water_table_elevation is not None and saturated_thickness is not None:
             raise RuntimeError(
                 "RadialUnconfinedDrawdown() must specify only "
-                + "water_table_elevation or saturated_thickness, but not "
-                + "both at the same time."
+                "water_table_elevation or saturated_thickness, but not "
+                "both at the same time."
             )
 
         if water_table_elevation is not None:
@@ -398,15 +398,14 @@ def _prop_check(self):
         if error:
             raise RuntimeError(
                 "RadialUnconfinedDrawdown: Attempted to solve radial "
-                + "groundwater model\nwith the following input not specified\n"
+                "groundwater model\nwith the following input not specified\n"
                 + "\n".join(error)
             )
         if self.well_top <= self.well_bot:
             raise RuntimeError(
                 "RadialUnconfinedDrawdown: "
-                + "well_screen_elevation_top <= well_screen_elevation_bottom\n"
-                + f"That is: {self.well_top} <= "
-                + f"{self.well_bot}"
+                "well_screen_elevation_top <= well_screen_elevation_bottom\n"
+                f"That is: {self.well_top} <= {self.well_bot}"
             )
 
     def drawdown(
@@ -581,7 +580,7 @@ def drawdown_times(
         if ty_time and ts_time:
             raise RuntimeError(
                 "RadialUnconfinedDrawdown.drawdown_times "
-                + "cannot set both ty_time and ts_time to True."
+                "cannot set both ty_time and ts_time to True."
             )
 
         r = radius
@@ -633,7 +632,7 @@ def drawdown_times(
         for stp, ts in enumerate(ts_list):
             if show_progress:
                 print(
-                    f"Solving {stp+1:4d} of {nstp}; " + f"time = {self.ts2time(ts, r)}",
+                    f"Solving {stp+1:4d} of {nstp}; time = {self.ts2time(ts, r)}",
                     end="",
                 )
 
@@ -721,23 +720,20 @@ def drawdown_times(
             bad_times = "\n".join([str(times[it]) for it in bessel_root_limit_reached])
             warnings.warn(
                 "\n\nRadialUnconfinedDrawdown.drawdown_times failed to "
-                + f"meet convergence sumrtol = {sumrtol}"
-                + "\nwithin the precalculated Bessel root solutions "
-                + "(convergence is evaluated at every second Bessel root).\n\n"
-                + "The number of Bessel roots are automatically increased "
-                + "up to:\n"
-                + f"   {bessel_roots0} * 2^bessel_loop_limit\nwhere:\n"
-                + "   bessel_loop_limit = {bessel_loop_limit}\n"
-                + f"resulting in {1024*2**bessel_loop_limit} roots evaluated, "
-                + "with the last root being {root}\n"
-                + f"(That is, the Neuman integral was solved form 0 to {root})"
-                + "\n\n"
-                + "You can either ignore this warning\n"
-                + "or to remove it attempt to increase bessel_loop_limit\n"
-                + "or increase sumrtol (reducing accuracy).\n\nThe following "
-                + "times are what triggered this warning:\n"
-                + bad_times
-                + "\n"
+                f"meet convergence sumrtol = {sumrtol}"
+                "\nwithin the precalculated Bessel root solutions "
+                "(convergence is evaluated at every second Bessel root).\n\n"
+                "The number of Bessel roots are automatically increased "
+                "up to:\n"
+                f"   {bessel_roots0} * 2^bessel_loop_limit\nwhere:\n"
+                "   bessel_loop_limit = {bessel_loop_limit}\n"
+                f"resulting in {1024*2**bessel_loop_limit} roots evaluated, "
+                "with the last root being {root}\n"
+                f"(That is, the Neuman integral was solved form 0 to {root})\n\n"
+                "You can either ignore this warning\n"
+                "or to remove it attempt to increase bessel_loop_limit\n"
+                "or increase sumrtol (reducing accuracy).\n\nThe following "
+                "times are what triggered this warning:\n" + bad_times + "\n"
             )
         return s
 
@@ -1012,7 +1008,7 @@ def _get_bracket(func, a, b, arg=(), internal_search_split=100):
 
         raise RuntimeError(
             "get_bracket: failed to find bracket interval with opposite "
-            + f"signs, that is: f(a)*f(b) < 0 for func: {func}"
+            f"signs, that is: f(a)*f(b) < 0 for func: {func}"
         )
 
 
@@ -1150,7 +1146,7 @@ def get_radius_lay_from_node(node, nradial):
     ["h_mid", "head", (get_radial_node(11, (nlay - 1) // 2, nradial),)],
     ["h_bot", "head", (get_radial_node(11, nlay - 1, nradial),)],
 ]
-obsdict = {"{}.obs.head.csv".format(sim_name): obslist}
+obsdict = {f"{sim_name}.obs.head.csv": obslist}
 
 # Solver parameters
 
@@ -1827,7 +1823,7 @@ def plot_ts(sim, verbose=False, solve_analytical_solution=False):
         fig.tight_layout()
 
         if plot_save:
-            fpth = figs_path / "{}-{}{}".format(sim_name, obs_fig, ".png")
+            fpth = figs_path / f"{sim_name}-{obs_fig}.png"
             fig.savefig(fpth)
 
         obs_fig = "obs-dimensionless"
@@ -1869,7 +1865,7 @@ def plot_ts(sim, verbose=False, solve_analytical_solution=False):
         fig.tight_layout()
 
         if plot_save:
-            fpth = figs_path / "{}-{}{}".format(sim_name, obs_fig, ".png")
+            fpth = figs_path / f"{sim_name}-{obs_fig}.png"
             fig.savefig(fpth)
 
 
@@ -1943,7 +1939,7 @@ def plot_grid(verbose=False):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}-grid{}".format(sim_name, ".png")
+            fpth = figs_path / f"{sim_name}-grid.png"
             fig.savefig(fpth)
 
 
diff --git a/scripts/ex-gwf-sagehen.py b/scripts/ex-gwf-sagehen.py
index 472ffa55..0ac4db01 100644
--- a/scripts/ex-gwf-sagehen.py
+++ b/scripts/ex-gwf-sagehen.py
@@ -2166,11 +2166,11 @@ def plot_results(mf6):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(mf6.name + "-finfFact", ".png")
+            fpth = figs_path / f"{mf6.name}-finfFact.png"
             fig.savefig(fpth)
 
         # Start by retrieving some output
-        gwf = mf6.get_model(list(mf6.model_names)[0])
+        gwf = mf6.get_model(next(iter(mf6.model_names)))
         hdsobj = gwf.output.head()
         modobj = gwf.output.budget()
         sfrobj = gwf.sfr.output.budget()
@@ -2204,7 +2204,7 @@ def plot_results(mf6):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(mf6.name + "-gwDepth", ".png")
+            fpth = figs_path / f"{mf6.name}-gwDepth.png"
             fig.savefig(fpth)
 
         drn_disQ = []
@@ -2332,7 +2332,7 @@ def plot_results(mf6):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(mf6.name + "-uzFlow", ".png")
+            fpth = figs_path / f"{mf6.name}-uzFlow.png"
             fig.savefig(fpth)
 
         data_sfr = {
@@ -2374,7 +2374,7 @@ def plot_results(mf6):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(mf6.name + "-swFlow", ".png")
+            fpth = figs_path / f"{mf6.name}-swFlow.png"
             fig.savefig(fpth)
 
 
diff --git a/scripts/ex-gwf-twri.py b/scripts/ex-gwf-twri.py
index 6df6c6ff..a336f572 100644
--- a/scripts/ex-gwf-twri.py
+++ b/scripts/ex-gwf-twri.py
@@ -300,10 +300,7 @@ def plot_results(sim, mf, silent=True):
         # check that the results are comparable
         for idx, k in enumerate((0, 2, 4)):
             diff = np.abs(head[k] - head0[idx])
-            msg = (
-                "aquifer {}: ".format(idx + 1)
-                + f"maximum absolute head difference is {diff.max()}"
-            )
+            msg = f"aquifer {idx + 1}: maximum absolute head difference is {diff.max()}"
             assert diff.max() < 0.05, msg
             if not silent:
                 print(msg)
diff --git a/scripts/ex-gwt-gwtgwt-p10.py b/scripts/ex-gwt-gwtgwt-p10.py
index 17274147..bd57dbe3 100644
--- a/scripts/ex-gwt-gwtgwt-p10.py
+++ b/scripts/ex-gwt-gwtgwt-p10.py
@@ -247,14 +247,14 @@ def build_models():
         exgtype="GWF6-GWT6",
         exgmnamea=gwfname_out,
         exgmnameb=gwtname_out,
-        filename="{}.gwfgwt".format("outer"),
+        filename="outer.gwfgwt",
     )
     flopy.mf6.ModflowGwfgwt(
         sim,
         exgtype="GWF6-GWT6",
         exgmnamea=gwfname_inn,
         exgmnameb=gwtname_inn,
-        filename="{}.gwfgwt".format("inner"),
+        filename="inner.gwfgwt",
     )
 
     return sim
@@ -277,7 +277,7 @@ def add_flow(sim):
         scaling_method="NONE",
         reordering_method="NONE",
         relaxation_factor=relax,
-        filename="{}.ims".format("gwfsolver"),
+        filename="gwfsolver.ims",
     )
 
     gwf_outer = add_outer_gwfmodel(sim)
@@ -576,7 +576,7 @@ def add_transport(sim):
         scaling_method="NONE",
         reordering_method="NONE",
         relaxation_factor=relax,
-        filename="{}.ims".format("gwtsolver"),
+        filename="gwtsolver.ims",
     )
 
     # Instantiating transport advection package
diff --git a/scripts/ex-gwt-hecht-mendez.py b/scripts/ex-gwt-hecht-mendez.py
index 2706b959..f5ff8ca5 100644
--- a/scripts/ex-gwt-hecht-mendez.py
+++ b/scripts/ex-gwt-hecht-mendez.py
@@ -938,7 +938,7 @@ def plot_results(
             plt.show()
         if plot_save:
             letter = chr(ord("@") + idx + 1)
-            fpth = figs_path / "{}{}".format("ex-" + sim_name + "-" + letter, ".png")
+            fpth = figs_path / f"ex-{sim_name}-{letter}.png"
             fig.savefig(fpth)
 
 
diff --git a/scripts/ex-gwt-keating.py b/scripts/ex-gwt-keating.py
index 116b01ff..d9163ea4 100644
--- a/scripts/ex-gwt-keating.py
+++ b/scripts/ex-gwt-keating.py
@@ -489,7 +489,7 @@ def update(i):
         idx_end = (np.abs(conc_times - 18000.0)).argmin()
         ani = FuncAnimation(fig, update, range(1, idx_end), init_func=init)
         writer = PillowWriter(fps=25)
-        fpth = figs_path / "{}{}".format(sim_name, ".gif")
+        fpth = figs_path / f"{sim_name}.gif"
         ani.save(fpth, writer=writer)
 
 
diff --git a/scripts/ex-gwt-mt3dms-p01.py b/scripts/ex-gwt-mt3dms-p01.py
index 31605617..498461bd 100644
--- a/scripts/ex-gwt-mt3dms-p01.py
+++ b/scripts/ex-gwt-mt3dms-p01.py
@@ -560,7 +560,7 @@ def plot_results(mt3d, mf6, idx, ax=None):
         ax.set_xlim(0, 1000)
         ax.set_xlabel("Distance, in m")
         ax.set_ylabel("Concentration")
-        title = "Concentration Profile at Time = 2,000 " + "{}".format(time_units)
+        title = f"Concentration Profile at Time = 2,000 {time_units}"
         ax.legend()
         letter = chr(ord("@") + idx + 1)
         styles.heading(letter=letter, heading=title)
diff --git a/scripts/ex-gwt-mt3dms-p03.py b/scripts/ex-gwt-mt3dms-p03.py
index 1a583bb3..06272954 100644
--- a/scripts/ex-gwt-mt3dms-p03.py
+++ b/scripts/ex-gwt-mt3dms-p03.py
@@ -518,7 +518,7 @@ def plot_results(mt3d, mf6, idx, ax=None):
 
         plt.xlabel("DISTANCE ALONG X-AXIS, IN METERS")
         plt.ylabel("DISTANCE ALONG Y-AXIS, IN METERS")
-        title = "Plume at Time = 365 " + f"{time_units}"
+        title = f"Plume at Time = 365 {time_units}"
 
         ax.legend(lines, labels, loc="upper left")
 
diff --git a/scripts/ex-gwt-mt3dms-p05.py b/scripts/ex-gwt-mt3dms-p05.py
index 6f9ed859..74683393 100644
--- a/scripts/ex-gwt-mt3dms-p05.py
+++ b/scripts/ex-gwt-mt3dms-p05.py
@@ -527,7 +527,7 @@ def plot_results(mt3d, mf6, idx, ax=None, ax2=None):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(sim_name + "-xsec", ".png")
+            fpth = figs_path / f"{sim_name}-xsec.png"
             fig.savefig(fpth)
 
         # second plot
@@ -565,7 +565,7 @@ def plot_results(mt3d, mf6, idx, ax=None, ax2=None):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(sim_name + "-planView", ".png")
+            fpth = figs_path / f"{sim_name}-planView.png"
             fig.savefig(fpth)
 
 
diff --git a/scripts/ex-gwt-mt3dms-p08.py b/scripts/ex-gwt-mt3dms-p08.py
index 38dcc99d..eed2d29e 100644
--- a/scripts/ex-gwt-mt3dms-p08.py
+++ b/scripts/ex-gwt-mt3dms-p08.py
@@ -627,9 +627,7 @@ def plot_results(mf2k5, mt3d, mf6, idx, ax=None):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(
-                mf6.name + "-" + str(yr_idx[i] + 1) + "yrs", ".png"
-            )
+            fpth = figs_path / f"{mf6.name}-{yr_idx[i] + 1}yrs.png"
             fig.savefig(fpth)
 
         # Plot after 12 years
@@ -671,9 +669,7 @@ def plot_results(mf2k5, mt3d, mf6, idx, ax=None):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(
-                mf6.name + "-" + str(yr_idx[i] + 1) + "yrs", ".png"
-            )
+            fpth = figs_path / f"{mf6.name}-{yr_idx[i] + 1}yrs.png"
             fig.savefig(fpth)
 
         # Plot after 20 years
@@ -715,9 +711,7 @@ def plot_results(mf2k5, mt3d, mf6, idx, ax=None):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(
-                mf6.name + "-" + str(yr_idx[i] + 1) + "yrs", ".png"
-            )
+            fpth = figs_path / f"{mf6.name}-{yr_idx[i] + 1}yrs.png"
             fig.savefig(fpth)
 
 
diff --git a/scripts/ex-gwt-mt3dms-p09.py b/scripts/ex-gwt-mt3dms-p09.py
index 5a155566..800a7441 100644
--- a/scripts/ex-gwt-mt3dms-p09.py
+++ b/scripts/ex-gwt-mt3dms-p09.py
@@ -579,7 +579,7 @@ def plot_results(mf2k5, mt3d, mf6, idx, ax=None):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(sim_name, ".png")
+            fpth = figs_path / f"{sim_name}.png"
             fig.savefig(fpth)
 
 
diff --git a/scripts/ex-gwt-mt3dms-p10.py b/scripts/ex-gwt-mt3dms-p10.py
index c9b5a812..06ef58a6 100644
--- a/scripts/ex-gwt-mt3dms-p10.py
+++ b/scripts/ex-gwt-mt3dms-p10.py
@@ -749,7 +749,7 @@ def plot_results(mf2k5, mt3d, mf6, idx, ax=None):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(mf6.name, ".png")
+            fpth = figs_path / f"{mf6.name}.png"
             fig.savefig(fpth)
 
 
diff --git a/scripts/ex-gwt-mt3dsupp632.py b/scripts/ex-gwt-mt3dsupp632.py
index 24da3d40..ac71f2d2 100644
--- a/scripts/ex-gwt-mt3dsupp632.py
+++ b/scripts/ex-gwt-mt3dsupp632.py
@@ -425,7 +425,7 @@ def plot_results():
         if plot_show:
             plt.show()
         if plot_save:
-            fname = "{}{}".format("ex-gwt-mt3dsupp632", ".png")
+            fname = "ex-gwt-mt3dsupp632.png"
             fpth = figs_path / fname
             fig.savefig(fpth)
 
diff --git a/scripts/ex-gwt-rotate.py b/scripts/ex-gwt-rotate.py
index 08dcc991..a618e0fe 100644
--- a/scripts/ex-gwt-rotate.py
+++ b/scripts/ex-gwt-rotate.py
@@ -479,7 +479,7 @@ def update(i):
 
         ani = FuncAnimation(fig, update, range(1, times.shape[0], 5), init_func=init)
         writer = PillowWriter(fps=50)
-        fpth = figs_path / "{}{}".format(sim_name, ".gif")
+        fpth = figs_path / f"{sim_name}.gif"
         ani.save(fpth, writer=writer)
 
 
diff --git a/scripts/ex-gwt-saltlake.py b/scripts/ex-gwt-saltlake.py
index 2ceaa0b3..7939f1fd 100644
--- a/scripts/ex-gwt-saltlake.py
+++ b/scripts/ex-gwt-saltlake.py
@@ -344,7 +344,7 @@ def update(i):
 
         ani = FuncAnimation(fig, update, range(1, times.shape[0]), init_func=init)
         writer = PillowWriter(fps=50)
-        fpth = figs_path / "{}{}".format(sim_name, ".gif")
+        fpth = figs_path / f"{sim_name}.gif"
         ani.save(fpth, writer=writer)
 
 
diff --git a/scripts/ex-gwt-stallman.py b/scripts/ex-gwt-stallman.py
index f905eaec..af5dcbe5 100644
--- a/scripts/ex-gwt-stallman.py
+++ b/scripts/ex-gwt-stallman.py
@@ -387,7 +387,7 @@ def update(j):
         return line
 
     ani = animation.FuncAnimation(fig, update, times.shape[0], init_func=init)
-    fpth = figs_path / "{}{}".format(sim_name, ".gif")
+    fpth = figs_path / f"{sim_name}.gif"
     ani.save(fpth, fps=50)
 
 
diff --git a/scripts/ex-gwt-uzt-2d.py b/scripts/ex-gwt-uzt-2d.py
index 90bbcd3a..bd95e591 100644
--- a/scripts/ex-gwt-uzt-2d.py
+++ b/scripts/ex-gwt-uzt-2d.py
@@ -927,7 +927,7 @@ def plot_results(mf2k5, mt3d, mf6, idx, ax=None):
         if plot_show:
             plt.show()
         if plot_save:
-            fpth = figs_path / "{}{}".format(sim_name, ".png")
+            fpth = figs_path / f"{sim_name}.png"
             fig.savefig(fpth)
 
 
diff --git a/scripts/ex-prt-mp7-p01.py b/scripts/ex-prt-mp7-p01.py
index e8649ac0..8caf01a5 100644
--- a/scripts/ex-prt-mp7-p01.py
+++ b/scripts/ex-prt-mp7-p01.py
@@ -622,7 +622,7 @@ def plot_head(gwf, head):
         cint = 0.25
         hmin = head[ilay, 0, :].min()
         hmax = head[ilay, 0, :].max()
-        styles.heading(ax=ax, heading=f"Head, layer {str(ilay + 1)}")
+        styles.heading(ax=ax, heading=f"Head, layer {ilay + 1!s}")
         mm = flopy.plot.PlotMapView(gwf, ax=ax, layer=ilay)
         mm.plot_grid(lw=0.5)
         mm.plot_bc("WEL", plotAll=True)
diff --git a/scripts/ex-prt-mp7-p02.py b/scripts/ex-prt-mp7-p02.py
index 9fc31507..041679e3 100644
--- a/scripts/ex-prt-mp7-p02.py
+++ b/scripts/ex-prt-mp7-p02.py
@@ -301,7 +301,7 @@ def build_gwf_sim():
 
     # Instantiate the MODFLOW 6 gwf (groundwater-flow) model
     gwf = flopy.mf6.ModflowGwf(
-        sim, modelname=gwf_name, model_nam_file="{}.nam".format(gwf_name)
+        sim, modelname=gwf_name, model_nam_file=f"{gwf_name}.nam"
     )
     gwf.name_file.save_flows = True
 
@@ -355,9 +355,9 @@ def build_gwf_sim():
     )
 
     # Instantiate the MODFLOW 6 gwf output control package
-    headfile = "{}.hds".format(gwf_name)
+    headfile = f"{gwf_name}.hds"
     head_record = [headfile]
-    budgetfile = "{}.cbb".format(gwf_name)
+    budgetfile = f"{gwf_name}.cbb"
     budget_record = [budgetfile]
     flopy.mf6.ModflowGwfoc(
         gwf,
@@ -406,7 +406,7 @@ def build_prt_sim():
 
     # Instantiate the MODFLOW 6 prt model
     prt = flopy.mf6.ModflowPrt(
-        sim, modelname=prt_name, model_nam_file="{}.nam".format(prt_name)
+        sim, modelname=prt_name, model_nam_file=f"{prt_name}.nam"
     )
 
     # Instantiate the MODFLOW 6 prt discretization package
@@ -443,9 +443,9 @@ def add_prp(part):
     add_prp("B")
 
     # Instantiate the MODFLOW 6 prt output control package
-    budgetfile = "{}.bud".format(prt_name)
-    trackfile = "{}.trk".format(prt_name)
-    trackcsvfile = "{}.trk.csv".format(prt_name)
+    budgetfile = f"{prt_name}.bud"
+    trackfile = f"{prt_name}.trk"
+    trackcsvfile = f"{prt_name}.trk.csv"
     budget_record = [budgetfile]
     track_record = [trackfile]
     trackcsv_record = [trackcsvfile]
@@ -473,7 +473,7 @@ def add_prp(part):
     ems = flopy.mf6.ModflowEms(
         sim,
         pname="ems",
-        filename="{}.ems".format(prt_name),
+        filename=f"{prt_name}.ems",
     )
     sim.register_solution_package(ems, [prt.name])
 
@@ -856,7 +856,7 @@ def plot_head(gwf, head, ibd):
         cint = 0.25
         hmin = head[ilay, 0, :].min()
         hmax = head[ilay, 0, :].max()
-        styles.heading(ax=ax, heading="Head, layer {}".format(ilay + 1))
+        styles.heading(ax=ax, heading=f"Head, layer {ilay + 1}")
         mm = flopy.plot.PlotMapView(gwf, ax=ax, layer=ilay)
         mm.plot_bc("WEL", plotAll=True)
         mm.plot_bc("RIV", plotAll=True, color="teal")
@@ -890,7 +890,7 @@ def plot_head(gwf, head, ibd):
         if plot_show:
             plt.show()
         if plot_save:
-            fig.savefig(figs_path / "{}-head".format(sim_name))
+            fig.savefig(figs_path / f"{sim_name}-head")
 
 
 def plot_points(ax, gwf, data):
@@ -992,7 +992,7 @@ def plot_pathlines_and_points(gwf, mf6pl, ibd, title=None):
     if plot_show:
         plt.show()
     if plot_save:
-        fig.savefig(figs_path / "{}-paths".format(sim_name))
+        fig.savefig(figs_path / f"{sim_name}-paths")
 
 
 def plot_2b(gwf, mf6endpoints, ibd, title=None):
@@ -1011,7 +1011,7 @@ def plot_2b(gwf, mf6endpoints, ibd, title=None):
     if plot_show:
         plt.show()
     if plot_save:
-        fig.savefig(figs_path / "{}-endpts".format(sim_name))
+        fig.savefig(figs_path / f"{sim_name}-endpts")
 
 
 def plot_3d(gwf, pathlines, endpoints=None, title=None):
diff --git a/scripts/ex-prt-mp7-p03.py b/scripts/ex-prt-mp7-p03.py
index f66b34ca..b71cd1e7 100644
--- a/scripts/ex-prt-mp7-p03.py
+++ b/scripts/ex-prt-mp7-p03.py
@@ -373,7 +373,7 @@ def build_prt_model():
 
     # Instantiate the MODFLOW 6 prt model
     prt = flopy.mf6.ModflowPrt(
-        sim, modelname=prt_name, model_nam_file="{}.nam".format(prt_name)
+        sim, modelname=prt_name, model_nam_file=f"{prt_name}.nam"
     )
 
     # Instantiate the MODFLOW 6 prt discretization package
@@ -404,7 +404,7 @@ def build_prt_model():
     flopy.mf6.ModflowPrtprp(
         prt,
         pname="prp1",
-        filename="{}_1.prp".format(prt_name),
+        filename=f"{prt_name}_1.prp",
         nreleasepts=len(release_points),
         packagedata=release_points,
         nreleasetimes=len(release_times),
@@ -641,7 +641,7 @@ def plot_head(gwf, head):
         cint = 0.25
         hmin = head[ilay, 0, :].min()
         hmax = head[ilay, 0, :].max()
-        styles.heading(ax=ax, heading=f"Head, layer {str(ilay + 1)}, time=0")
+        styles.heading(ax=ax, heading=f"Head, layer {ilay + 1!s}, time=0")
         mm = flopy.plot.PlotMapView(gwf, ax=ax, layer=ilay)
         mm.plot_grid(lw=0.5)
         mm.plot_bc("WEL", plotAll=True)
diff --git a/scripts/ex-prt-mp7-p04.py b/scripts/ex-prt-mp7-p04.py
index 7c3d290f..b8fb67e0 100644
--- a/scripts/ex-prt-mp7-p04.py
+++ b/scripts/ex-prt-mp7-p04.py
@@ -570,7 +570,7 @@ def build_prt():
 
     # Instantiate the MODFLOW 6 prt model
     prt = flopy.mf6.ModflowPrt(
-        simprt, modelname=prt_name, model_nam_file="{}.nam".format(prt_name)
+        simprt, modelname=prt_name, model_nam_file=f"{prt_name}.nam"
     )
 
     # Instantiate the MODFLOW 6 DISV vertex grid discretization
@@ -583,7 +583,7 @@ def build_prt():
     flopy.mf6.ModflowPrtprp(
         prt,
         pname="prp",
-        filename="{}_4.prp".format(prt_name),
+        filename=f"{prt_name}_4.prp",
         nreleasepts=len(particles_prt),
         packagedata=particles_prt,
         perioddata={0: ["FIRST"]},
@@ -592,10 +592,10 @@ def build_prt():
     )
 
     # Instantiate the MODFLOW 6 prt output control package
-    budgetfile_prt = "{}.cbb".format(prt_name)
+    budgetfile_prt = f"{prt_name}.cbb"
     budget_record = [budgetfile_prt]
-    trackfile_prt = "{}.trk".format(prt_name)
-    trackcsvfile_prt = "{}.trk.csv".format(prt_name)
+    trackfile_prt = f"{prt_name}.trk"
+    trackcsvfile_prt = f"{prt_name}.trk.csv"
     flopy.mf6.ModflowPrtoc(
         prt,
         pname="oc",
@@ -618,7 +618,7 @@ def build_prt():
     ems = flopy.mf6.ModflowEms(
         simprt,
         pname="ems",
-        filename="{}.ems".format(prt_name),
+        filename=f"{prt_name}.ems",
     )
     simprt.register_solution_package(ems, [prt.name])
 
diff --git a/scripts/process-scripts.py b/scripts/process-scripts.py
index 93bc5e27..3c2c56b6 100644
--- a/scripts/process-scripts.py
+++ b/scripts/process-scripts.py
@@ -159,11 +159,9 @@ def make_tables():
                         except:
                             units = " (unknown)"
                     if len(table_line) > 0:
-                        table_line += "\t{}{} & {}".format(text_to_write, units, value)
+                        table_line += f"\t{text_to_write}{units} & {value}"
                     else:
-                        table_line = "\t& & {}{} & {}".format(
-                            text_to_write, units, value
-                        )
+                        table_line = f"\t& & {text_to_write}{units} & {value}"
                     table_line += " \\\\\n"
                     if len(row_color) > 0:
                         f.write(row_color)
@@ -369,7 +367,7 @@ def build_md_tables(ex_dict):
     for ex_name in ex_paks.keys():
         for ex_root in ex_order:
             if ex_root in ex_name:
-                pak_link[ex_name] = "[{}](_examples/{}.html)".format(ex_name, ex_root)
+                pak_link[ex_name] = f"[{ex_name}](_examples/{ex_root}.html)"
                 break
         if ex_name not in list(pak_link.keys()):
             pak_link[ex_name] = ex_name
@@ -396,23 +394,21 @@ def build_md_tables(ex_dict):
         line = "### Introduction\n\n"
         line += (
             "This document describes example problems for MODFLOW 6. The examples "
-            + "demonstrate use of the capabilities for select components of "
-            + "MODFLOW 6. A pdf of the examples can be downloaded from "
-            + "[ReadtheDocs](https://modflow6-examples.readthedocs.io/en/latest/) "
-            + "or from the [current release](https://github.com/MODFLOW-USGS/"
-            + "modflow6-examples/releases/download/current/mf6examples.pdf) on "
-            + "[GitHub](https://github.com/MODFLOW-USGS/modflow6-examples/).\n\n"
-            + "Examples have been included for the MODFLOW 6 components "
-            + "summarized in the tables below.\n\n"
+            "demonstrate use of the capabilities for select components of "
+            "MODFLOW 6. A pdf of the examples can be downloaded from "
+            "[ReadtheDocs](https://modflow6-examples.readthedocs.io/en/latest/) "
+            "or from the [current release](https://github.com/MODFLOW-USGS/"
+            "modflow6-examples/releases/download/current/mf6examples.pdf) on "
+            "[GitHub](https://github.com/MODFLOW-USGS/modflow6-examples/).\n\n"
+            "Examples have been included for the MODFLOW 6 components "
+            "summarized in the tables below.\n\n"
         )
         f.write(line)
 
         join_fmt = ", "
         header = (
-            "| Package | Examples                                            "
-            + "               |\n"
-            + "|---------|------------------------------"
-            + "--------------------------------------|\n"
+            "| Package | Examples                                                           |\n"
+            "|---------|--------------------------------------------------------------------|\n"
         )
         footer = "\n\n"