diff --git a/exercises/gwf_advanced/gwf-uzf.ipynb b/exercises/gwf_advanced/gwf-uzf.ipynb
new file mode 100644
index 0000000..6d80729
--- /dev/null
+++ b/exercises/gwf_advanced/gwf-uzf.ipynb
@@ -0,0 +1,821 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "96e0c173-fa09-49dd-9550-9e1c40d481b7",
+   "metadata": {},
+   "source": [
+    "## Exploring UZF\n",
+    "\n",
+    "This worked example sets up a MODFLOW 6 simulation with multiple GWF models to explore the sensitivity of various UZF parameters (e.g., Brooks-Corey epsilon, vks, etc.) within the UZF package and what impact it can have on the flow solution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "814a0872-52da-4c23-aa47-137635cc291e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Imports\n",
+    "import os\n",
+    "import sys\n",
+    "\n",
+    "# Suppress warning messages\n",
+    "import warnings\n",
+    "from pathlib import Path\n",
+    "\n",
+    "import flopy\n",
+    "import matplotlib\n",
+    "import matplotlib.animation as animation\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "from IPython.display import HTML\n",
+    "from matplotlib.collections import LineCollection\n",
+    "\n",
+    "warnings.simplefilter(\"ignore\", UserWarning)\n",
+    "warnings.simplefilter(\"ignore\", DeprecationWarning)\n",
+    "warnings.simplefilter(\"ignore\", FutureWarning)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "11ecceab-b4d2-4fcd-940e-cf9178f45f75",
+   "metadata": {},
+   "source": [
+    "### Set simulation name and workspace"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f00ee1e8-b907-40b6-83fe-98e8a41d3e3b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim_name = \"unsat-comp\"\n",
+    "sim_ws = Path(\"./temp/gwf-uzf\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6517e86c-a122-494d-9e18-a090c2c91213",
+   "metadata": {},
+   "source": [
+    "### Set some model input values"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cf334b1a-b9be-454d-93bb-a36be8204c25",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Model units\n",
+    "length_units = \"meters\"\n",
+    "time_units = \"seconds\"\n",
+    "\n",
+    "# Solver parameters\n",
+    "nouter, ninner = 100, 300\n",
+    "hclose, rclose, relax = 1e-8, 1e-8, 0.97\n",
+    "\n",
+    "# Set some parameter values\n",
+    "nper = 365  # Number of periods\n",
+    "nstp = 1  # Number of time steps\n",
+    "perlen = 86400  # Simulation time length ($s$)\n",
+    "nlay = 120  # Number of layers\n",
+    "nrow = 1  # Number of rows\n",
+    "ncol = 1  # Number of columns\n",
+    "system_length = 60.0  # Length of system ($m$)\n",
+    "delr = 1.0  # Column width ($m$)\n",
+    "delc = 1.0  # Row width ($m$)\n",
+    "delv_str = \"ranges from 0.1 to 1\"  # Layer thickness\n",
+    "top = 60.0  # Top of the model ($m$)\n",
+    "keff = 3e-6\n",
+    "hydraulic_conductivity = keff  # Hydraulic conductivity ($m s^{-1}$)\n",
+    "\n",
+    "Zw = 55\n",
+    "\n",
+    "thtr = 0.01\n",
+    "n1 = 0.24\n",
+    "n2 = 0.14"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ca534969-e061-408e-b834-ebf519a8a938",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Stress period input\n",
+    "per_data = []\n",
+    "for k in range(nper):\n",
+    "    per_data.append((perlen, nstp, 1.0))\n",
+    "per_mf6 = per_data\n",
+    "\n",
+    "# Geometry input\n",
+    "tp = top\n",
+    "botm = []\n",
+    "for i in range(nlay):\n",
+    "    if i == 0:\n",
+    "        botm.append(59.9)\n",
+    "    elif i == 119:\n",
+    "        botm.append(0.0)\n",
+    "    else:\n",
+    "        botm.append(60 - i * 0.5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "91398dd4-8a35-4f2f-9a98-c09a14994fbc",
+   "metadata": {},
+   "source": [
+    "### Setting some boundary package data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0e41a4d9-452f-48ac-add1-626f6b8dbfd2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Head input\n",
+    "chd_data = {}\n",
+    "chd_data[0] = [[(119, 0, 0), 53.875]]\n",
+    "\n",
+    "chd_mf6 = chd_data\n",
+    "\n",
+    "# Recharge input\n",
+    "wel_data = {}\n",
+    "wel_data[0] = [[(0, 0, 0), \"prcp\"]]\n",
+    "wel_mf6 = wel_data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5c75f711-a5e2-4e01-9661-9b796c3096a9",
+   "metadata": {},
+   "source": [
+    "### Read in time series input and prepare it for MODFLOW 6 time series utility"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "340e1f3d-7f03-45b3-979c-97543367f06d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Boundary temperature input\n",
+    "data_path = \"../../data/uzf\"\n",
+    "df = pd.read_csv(os.path.join(data_path, \"Daily_dat.csv\"), index_col=0)\n",
+    "\n",
+    "ctp_spd = {}\n",
+    "for i in range(nper):\n",
+    "    ctp_spd[i] = [\n",
+    "        [(0, 0, 0), df[\"TAVG_rmFrz\"][i]]\n",
+    "    ]  # ctp_spd.update({i: [[(0, 0, 0), df['TAVG_rmFrz'][i]]]})\n",
+    "\n",
+    "prcp_ts_data = []\n",
+    "pet_ts_data = []\n",
+    "temp_ts_data = []\n",
+    "\n",
+    "for n in range(0, len(df)):\n",
+    "    time = n * 86400\n",
+    "    prcp = df[\"PRCP\"][n] / 1000 / 86400\n",
+    "    pet = df[\"PET\"][n] / 1000 / 86400\n",
+    "    temp = df[\"TAVG_rmFrz\"][n]\n",
+    "    prcp_ts_data.append((time, prcp))\n",
+    "    temp_ts_data.append((time, temp))\n",
+    "    pet_ts_data.append((time, pet))\n",
+    "\n",
+    "prcp_ts_dict_base = {\n",
+    "    \"filename\": \"prcp.ts\",\n",
+    "    \"time_series_namerecord\": \"prcp\",\n",
+    "    \"timeseries\": prcp_ts_data,\n",
+    "    \"interpolation_methodrecord\": \"linearend\",\n",
+    "}\n",
+    "\n",
+    "# Make another copy of the time series file since two GWF model's cannot use the same input file (a file unit number issue)\n",
+    "prcp_ts_dict_eps = {\n",
+    "    \"filename\": \"prcp_eps.ts\",\n",
+    "    \"time_series_namerecord\": \"prcp_eps\",\n",
+    "    \"timeseries\": prcp_ts_data.copy(),\n",
+    "    \"interpolation_methodrecord\": \"linearend\",\n",
+    "}\n",
+    "\n",
+    "prcp_ts_dict_vks = {\n",
+    "    \"filename\": \"prcp_vks.ts\",\n",
+    "    \"time_series_namerecord\": \"prcp_vks\",\n",
+    "    \"timeseries\": prcp_ts_data.copy(),\n",
+    "    \"interpolation_methodrecord\": \"linearend\",\n",
+    "}\n",
+    "\n",
+    "prcp_ts_dict_uzet = {\n",
+    "    \"filename\": \"prcp_uzet.ts\",\n",
+    "    \"time_series_namerecord\": \"prcp_uzet\",\n",
+    "    \"timeseries\": prcp_ts_data.copy(),\n",
+    "    \"interpolation_methodrecord\": \"linearend\",\n",
+    "}\n",
+    "\n",
+    "prcp_ts_dict_theta_s = {\n",
+    "    \"filename\": \"prcp_theta_s.ts\",\n",
+    "    \"time_series_namerecord\": \"prcp_theta_s\",\n",
+    "    \"timeseries\": prcp_ts_data.copy(),\n",
+    "    \"interpolation_methodrecord\": \"linearend\",\n",
+    "}\n",
+    "\n",
+    "pet_ts_dict = {\n",
+    "    \"filename\": \"pet.ts\",\n",
+    "    \"time_series_namerecord\": \"pet\",\n",
+    "    \"timeseries\": pet_ts_data,\n",
+    "    \"interpolation_methodrecord\": \"linearend\",\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f6a14d50-fd64-42cb-90ba-fd60b40311be",
+   "metadata": {},
+   "source": [
+    "### Setup objects for UZF \n",
+    "\n",
+    "Includes two functions to help create the UZF input later on"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8eb8c6b5-781d-4cf2-bc7a-cedf10543e6e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# transient uzf info\n",
+    "def set_uzf_spd_info(gwfname, ietflg, surfdep):\n",
+    "    # The following is needed because different models within the same simulation cannot all access the same file.\n",
+    "    if \"base\" in gwfname:\n",
+    "        finf = \"prcp\"\n",
+    "    elif \"eps\" in gwfname:\n",
+    "        finf = \"prcp_eps\"\n",
+    "    elif \"vks\" in gwfname:\n",
+    "        finf = \"prcp_vks\"\n",
+    "    elif \"uzet\" in gwfname:\n",
+    "        finf = \"prcp_uzet\"\n",
+    "    elif \"theta\" in gwfname:\n",
+    "        finf = \"prcp_theta_s\"\n",
+    "\n",
+    "    zero = 0.0\n",
+    "    extwc = 0.0\n",
+    "    if ietflg == 0:\n",
+    "        pet = 0.0\n",
+    "        extdp = 0.0\n",
+    "        uzf_spd = [[0, finf, pet, extdp, extwc, zero, zero, zero]]\n",
+    "    elif ietflg != 0:\n",
+    "        pet = \"pet\"\n",
+    "        extdp = 3.9  # meters\n",
+    "        uzf_spd = []\n",
+    "        for i in np.arange(nlay):\n",
+    "            if (top - surfdep) - botm[i] < extdp:\n",
+    "                if i > 0:\n",
+    "                    finf = 0.0\n",
+    "                uzf_spd.append([i, finf, pet, extdp, extwc, zero, zero, zero])\n",
+    "\n",
+    "    uzf_spd = {0: uzf_spd}\n",
+    "\n",
+    "    return uzf_spd\n",
+    "\n",
+    "\n",
+    "def add_uzf_obs(uz_obs_loc, obsnm, iuz, elev):\n",
+    "    uz_obs_loc.append((obsnm, \"water-content\", iuz, elev))\n",
+    "    return uz_obs_loc\n",
+    "\n",
+    "\n",
+    "def add_uzf_wc_profile_obs():\n",
+    "    uz_obs_loc = []\n",
+    "    iuz = 1\n",
+    "    depth = round((top - botm[0]) / 2, 2)\n",
+    "    elev = top - depth\n",
+    "    # Top uzf object (has a unique thickness)\n",
+    "    obsnm = \"uzf\" + str(iuz).zfill(2) + \"_elev=\" + str(elev)\n",
+    "    uz_obs_loc = add_uzf_obs(uz_obs_loc, obsnm, iuz, depth)\n",
+    "\n",
+    "    # Loop through upper X-number cells (17 cells currently)\n",
+    "    increment = 0.1\n",
+    "    for iuz in np.arange(1, 17, 1):\n",
+    "        depth1 = depth  # 0.05\n",
+    "        depth2 = round(depth1 + increment, 2)  # 0.15\n",
+    "        depth3 = round(depth2 + increment, 2)  # 0.25\n",
+    "        depth4 = round(depth3 + increment, 2)  # 0.35\n",
+    "        depth5 = round(depth4 + increment, 2)  # 0.45\n",
+    "        elev1 = round(botm[iuz - 1] - depth1, 2)\n",
+    "        elev2 = round(botm[iuz - 1] - depth2, 2)\n",
+    "        elev3 = round(botm[iuz - 1] - depth3, 2)\n",
+    "        elev4 = round(botm[iuz - 1] - depth4, 2)\n",
+    "        elev5 = round(botm[iuz - 1] - depth5, 2)\n",
+    "        obsnm1 = \"uzf\" + str(iuz).zfill(2) + \"_elev=\" + str(elev1)\n",
+    "        obsnm2 = \"uzf\" + str(iuz).zfill(2) + \"_elev=\" + str(elev2)\n",
+    "        obsnm3 = \"uzf\" + str(iuz).zfill(2) + \"_elev=\" + str(elev3)\n",
+    "        obsnm4 = \"uzf\" + str(iuz).zfill(2) + \"_elev=\" + str(elev4)\n",
+    "        obsnm5 = \"uzf\" + str(iuz).zfill(2) + \"_elev=\" + str(elev5)\n",
+    "        uz_obs_loc = add_uzf_obs(uz_obs_loc, obsnm1, iuz + 1, depth1)\n",
+    "        uz_obs_loc = add_uzf_obs(uz_obs_loc, obsnm2, iuz + 1, depth2)\n",
+    "        uz_obs_loc = add_uzf_obs(uz_obs_loc, obsnm3, iuz + 1, depth3)\n",
+    "        uz_obs_loc = add_uzf_obs(uz_obs_loc, obsnm4, iuz + 1, depth4)\n",
+    "        if iuz > 1:\n",
+    "            uz_obs_loc = add_uzf_obs(uz_obs_loc, obsnm5, iuz + 1, depth5)\n",
+    "\n",
+    "    return uz_obs_loc"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f6c4e8db-71ed-45ec-9a14-a6788902b129",
+   "metadata": {},
+   "source": [
+    "### Function for creating a GWF model\n",
+    "\n",
+    "Since multiple GWF models will be added to the MODFLOW 6 simulation, define a function to do the work multiple times.\n",
+    "\n",
+    "**Note**: The function call accepts arguments that define parameter values for the UZF package, for example the Brooks-Corey epsilon ($\\epsilon$). This allows different calls to the this function to instantiate the UZF package with different parameter specifications for exploring their impact on flow in the unsaturated zone."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e12d4e4d-4a50-47f8-9f9f-862c631e3d3e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def add_gwf_model(sim, gwfname, eps=4, keff=3e-6, ietflg=0):\n",
+    "    msg0 = \"EPSILON must be between 3.5 and 14.0\"\n",
+    "    assert eps > 3.5 and eps < 15.0, msg0\n",
+    "\n",
+    "    # Instantiating MODFLOW 6 groundwater flow model\n",
+    "    gwf = flopy.mf6.ModflowGwf(\n",
+    "        sim,\n",
+    "        modelname=gwfname,\n",
+    "        save_flows=True,\n",
+    "        newtonoptions=True,\n",
+    "        model_nam_file=\"{}.nam\".format(gwfname),\n",
+    "    )\n",
+    "\n",
+    "    # Instantiating MODFLOW 6 solver for flow model\n",
+    "    imsgwf = flopy.mf6.ModflowIms(\n",
+    "        sim,\n",
+    "        print_option=\"SUMMARY\",\n",
+    "        outer_dvclose=hclose,\n",
+    "        outer_maximum=nouter,\n",
+    "        under_relaxation=\"NONE\",\n",
+    "        inner_maximum=ninner,\n",
+    "        inner_dvclose=hclose,\n",
+    "        rcloserecord=rclose,\n",
+    "        linear_acceleration=\"BICGSTAB\",\n",
+    "        scaling_method=\"NONE\",\n",
+    "        reordering_method=\"NONE\",\n",
+    "        relaxation_factor=relax,\n",
+    "        filename=\"{}.ims\".format(gwfname),\n",
+    "    )\n",
+    "    sim.register_ims_package(imsgwf, [gwfname])\n",
+    "\n",
+    "    # Instantiating MODFLOW 6 discretization package\n",
+    "    flopy.mf6.ModflowGwfdis(\n",
+    "        gwf,\n",
+    "        length_units=length_units,\n",
+    "        nlay=nlay,\n",
+    "        nrow=nrow,\n",
+    "        ncol=ncol,\n",
+    "        delr=delr,\n",
+    "        delc=delc,\n",
+    "        top=top,\n",
+    "        botm=botm,\n",
+    "        filename=\"{}.dis\".format(gwfname),\n",
+    "    )\n",
+    "\n",
+    "    # Instantiating MODFLOW 6 node-property flow package\n",
+    "    flopy.mf6.ModflowGwfnpf(\n",
+    "        gwf,\n",
+    "        save_specific_discharge=True,\n",
+    "        icelltype=1,\n",
+    "        k=hydraulic_conductivity,\n",
+    "        filename=\"{}.npf\".format(gwfname),\n",
+    "    )\n",
+    "\n",
+    "    # Instantiate MODFLOW 6 storage package\n",
+    "    if \"theta\" in gwfname:\n",
+    "        n_use = n2\n",
+    "    else:\n",
+    "        n_use = n1\n",
+    "\n",
+    "    flopy.mf6.ModflowGwfsto(\n",
+    "        gwf,\n",
+    "        ss=1e-6,\n",
+    "        sy=n_use,\n",
+    "        iconvert=1,\n",
+    "        # steady_state={0: True},\n",
+    "        filename=\"{}.sto\".format(gwfname),\n",
+    "    )\n",
+    "\n",
+    "    # Instantiating MODFLOW 6 initial conditions package for flow model\n",
+    "    flopy.mf6.ModflowGwfic(gwf, strt=Zw, filename=\"{}.ic\".format(gwfname))\n",
+    "\n",
+    "    # Instantiating MODFLOW 6 constant head package\n",
+    "    flopy.mf6.ModflowGwfchd(\n",
+    "        gwf,\n",
+    "        stress_period_data=chd_mf6,\n",
+    "        filename=\"{}.chd\".format(gwfname),\n",
+    "    )\n",
+    "\n",
+    "    # Unsaturated-zone flow package\n",
+    "    uz_obs_loc = add_uzf_wc_profile_obs()\n",
+    "    uzf_obs = {f\"{gwfname}.uzfobs\": uz_obs_loc}\n",
+    "    if \"base\" in gwfname:\n",
+    "        prcp_ts_dict_entry = prcp_ts_dict_base\n",
+    "    elif \"eps\" in gwfname:\n",
+    "        prcp_ts_dict_entry = prcp_ts_dict_eps\n",
+    "    elif \"vks\" in gwfname:\n",
+    "        prcp_ts_dict_entry = prcp_ts_dict_vks\n",
+    "    elif \"uzet\" in gwfname:\n",
+    "        prcp_ts_dict_entry = prcp_ts_dict_uzet\n",
+    "    elif \"theta\" in gwfname:\n",
+    "        prcp_ts_dict_entry = prcp_ts_dict_theta_s\n",
+    "\n",
+    "    # The following is for the upper-most (land surface) cell\n",
+    "    # iuzno  cellid landflg ivertcn surfdp vks thtr thts thti eps [bndnm]\n",
+    "    surfdep = 0.00001\n",
+    "    uzf_pkdat = [\n",
+    "        [0, (0, 0, 0), 1, 1, surfdep, keff, thtr, thtr + n1, 0.025, eps]\n",
+    "    ]\n",
+    "\n",
+    "    # Continue building the UZF stack of objects\n",
+    "    for iuzno in np.arange(1, nlay, 1):\n",
+    "        if iuzno < nlay - 1:\n",
+    "            ivertconn = iuzno + 1\n",
+    "        else:\n",
+    "            ivertconn = -1\n",
+    "        # Add the UZF object to the packagedata list object\n",
+    "        uzf_pkdat.append(\n",
+    "            [\n",
+    "                iuzno,\n",
+    "                (iuzno, 0, 0),\n",
+    "                0,\n",
+    "                ivertconn,\n",
+    "                0.0,\n",
+    "                keff,\n",
+    "                thtr,\n",
+    "                thtr + n_use,\n",
+    "                0.025,\n",
+    "                eps,\n",
+    "            ]\n",
+    "        )\n",
+    "\n",
+    "    # Call function to return stress period ('PERIODDATA' block) information\n",
+    "    uzf_spd = set_uzf_spd_info(gwfname, ietflg, surfdep)\n",
+    "    if ietflg == 0:\n",
+    "        simulate_et = False\n",
+    "        unsat_etwc = False\n",
+    "        linear_gwet = False\n",
+    "    else:\n",
+    "        simulate_et = True\n",
+    "        unsat_etwc = True\n",
+    "        linear_gwet = True\n",
+    "\n",
+    "    uzf = flopy.mf6.ModflowGwfuzf(\n",
+    "        gwf,\n",
+    "        print_flows=True,\n",
+    "        save_flows=True,\n",
+    "        wc_filerecord=gwfname + \".uzfwc.bin\",\n",
+    "        timeseries=prcp_ts_dict_entry,\n",
+    "        simulate_et=simulate_et,\n",
+    "        unsat_etwc=unsat_etwc,\n",
+    "        simulate_gwseep=False,\n",
+    "        linear_gwet=linear_gwet,\n",
+    "        boundnames=False,\n",
+    "        observations=uzf_obs,\n",
+    "        ntrailwaves=15,\n",
+    "        nwavesets=500,\n",
+    "        nuzfcells=len(uzf_pkdat),\n",
+    "        packagedata=uzf_pkdat,\n",
+    "        perioddata=uzf_spd,\n",
+    "        budget_filerecord=f\"{gwfname}.uzf.bud\",\n",
+    "        pname=\"UZF-1\",\n",
+    "        filename=f\"{gwfname}.uzf\",\n",
+    "    )\n",
+    "    if \"uzet\" in gwfname:\n",
+    "        uzf.ts.append_package(\n",
+    "            filename=\"pet.ts\",\n",
+    "            time_series_namerecord=\"pet\",\n",
+    "            timeseries=pet_ts_data,\n",
+    "            interpolation_methodrecord=\"linearend\",\n",
+    "        )\n",
+    "\n",
+    "    # Instantiating MODFLOW 6 output control package for flow model\n",
+    "    flopy.mf6.ModflowGwfoc(\n",
+    "        gwf,\n",
+    "        head_filerecord=\"{}.hds\".format(gwfname),\n",
+    "        budget_filerecord=\"{}.cbc\".format(gwfname),\n",
+    "        headprintrecord=[(\"COLUMNS\", 10, \"WIDTH\", 15, \"DIGITS\", 6, \"GENERAL\")],\n",
+    "        saverecord=[(\"HEAD\", \"LAST\"), (\"BUDGET\", \"LAST\")],\n",
+    "        printrecord=[(\"HEAD\", \"LAST\"), (\"BUDGET\", \"LAST\")],\n",
+    "    )\n",
+    "\n",
+    "    return sim"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9a449816-d34c-490d-946f-717f4fa41eb5",
+   "metadata": {},
+   "source": [
+    "### Create the simulation and add the models"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "90a6f870-c8f0-4e01-8d1e-9f88c7f1026b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim = flopy.mf6.MFSimulation(\n",
+    "    sim_name=sim_name, sim_ws=sim_ws, exe_name=\"mf6\", version=\"mf6\"\n",
+    ")\n",
+    "\n",
+    "# Instantiating MODFLOW 6 time discretization\n",
+    "flopy.mf6.ModflowTdis(\n",
+    "    sim, nper=nper, perioddata=per_mf6, time_units=time_units\n",
+    ")\n",
+    "\n",
+    "# Add the first GWF model: \"base case\"\n",
+    "gwfnamebase = \"gwf-base\"\n",
+    "sim = add_gwf_model(sim, gwfnamebase, eps=3.6, keff=3e-6, ietflg=0)\n",
+    "\n",
+    "# Add a second GWF model: adjust Brooks-Corey epsilon\n",
+    "gwfnamescen1 = \"gwf-eps\"\n",
+    "sim = add_gwf_model(sim, gwfnamescen1, eps=7.2, keff=3e-6, ietflg=0)\n",
+    "\n",
+    "# Add a third GWF model: restore Brooks-Corey epsilon, adjust vertical hydraulic conductivity\n",
+    "gwfnamescen2 = \"gwf-vks\"\n",
+    "sim = add_gwf_model(\n",
+    "    sim, gwfnamescen2, eps=3.6, keff=3e-5, ietflg=0\n",
+    ")  # More like a clean sand\n",
+    "\n",
+    "# Add a fourth GWF model: params same as \"base case\" scenario but with ET\n",
+    "gwfnamescen3 = \"gwf-uzet\"\n",
+    "sim = add_gwf_model(sim, gwfnamescen3, ietflg=1)  # More like a clean sand\n",
+    "\n",
+    "# Add a fifth GWF model: explore effect of saturated water content\n",
+    "gwfnamescen4 = \"gwf-theta_s\"\n",
+    "sim = add_gwf_model(\n",
+    "    sim,\n",
+    "    gwfnamescen4,\n",
+    "    ietflg=0,\n",
+    ")  # More like a clean sand\n",
+    "\n",
+    "sim.write_simulation(silent=False)\n",
+    "\n",
+    "success, buff = sim.run_simulation(silent=False, report=True)\n",
+    "assert success, buff"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "67b38632-1664-44d8-8bb4-4c0710bab575",
+   "metadata": {},
+   "source": [
+    "### Load and animate UZF & UZE output"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c638f809-aaca-44a7-93e8-fb22c75702f1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# load temperature, model 1\n",
+    "gwf1 = sim.get_model(gwfnamebase)\n",
+    "gwf2 = sim.get_model(gwfnamescen1)\n",
+    "gwf3 = sim.get_model(gwfnamescen2)\n",
+    "gwf4 = sim.get_model(gwfnamescen3)\n",
+    "gwf5 = sim.get_model(gwfnamescen4)\n",
+    "\n",
+    "# Water content output\n",
+    "pth = sim.sim_path\n",
+    "flb = gwfnamebase + \".uzfobs\"\n",
+    "uzdat_base = pd.read_csv(os.path.join(pth, flb))\n",
+    "\n",
+    "fleps = gwfnamescen1 + \".uzfobs\"\n",
+    "uzdat_eps = pd.read_csv(os.path.join(pth, fleps))\n",
+    "\n",
+    "flvks = gwfnamescen2 + \".uzfobs\"\n",
+    "uzdat_vks = pd.read_csv(os.path.join(pth, flvks))\n",
+    "\n",
+    "fluzet = gwfnamescen3 + \".uzfobs\"\n",
+    "uzdat_uzet = pd.read_csv(os.path.join(pth, fluzet))\n",
+    "\n",
+    "fltheta = gwfnamescen4 + \".uzfobs\"\n",
+    "uzdat_theta = pd.read_csv(os.path.join(pth, fltheta))\n",
+    "\n",
+    "# Timing information for printing on animation\n",
+    "uzdat_base[\"time\"] = uzdat_base[\"time\"] / 86400\n",
+    "uzdat_base[\"time\"] = uzdat_base[\"time\"].astype(int)\n",
+    "uzdat_base.set_index(\"time\", inplace=True)\n",
+    "\n",
+    "uzdat_eps[\"time\"] = uzdat_eps[\"time\"] / 86400\n",
+    "uzdat_eps[\"time\"] = uzdat_eps[\"time\"].astype(int)\n",
+    "uzdat_eps.set_index(\"time\", inplace=True)\n",
+    "\n",
+    "uzdat_vks[\"time\"] = uzdat_vks[\"time\"] / 86400\n",
+    "uzdat_vks[\"time\"] = uzdat_vks[\"time\"].astype(int)\n",
+    "uzdat_vks.set_index(\"time\", inplace=True)\n",
+    "\n",
+    "uzdat_uzet[\"time\"] = uzdat_uzet[\"time\"] / 86400\n",
+    "uzdat_uzet[\"time\"] = uzdat_uzet[\"time\"].astype(int)\n",
+    "uzdat_uzet.set_index(\"time\", inplace=True)\n",
+    "\n",
+    "uzdat_theta[\"time\"] = uzdat_theta[\"time\"] / 86400\n",
+    "uzdat_theta[\"time\"] = uzdat_theta[\"time\"].astype(int)\n",
+    "uzdat_theta.set_index(\"time\", inplace=True)\n",
+    "\n",
+    "# Mine all observations from obs names\n",
+    "depths_base = uzdat_base.columns\n",
+    "depths_base = [\n",
+    "    round(float(d.strip().split(\"=\")[-1]), 2) for d in depths_base if \"=\" in d\n",
+    "]\n",
+    "\n",
+    "depths_eps = uzdat_eps.columns\n",
+    "depths_eps = [\n",
+    "    round(float(d.strip().split(\"=\")[-1]), 2) for d in depths_eps if \"=\" in d\n",
+    "]\n",
+    "\n",
+    "depths_vks = uzdat_vks.columns\n",
+    "depths_vks = [\n",
+    "    round(float(d.strip().split(\"=\")[-1]), 2) for d in depths_vks if \"=\" in d\n",
+    "]\n",
+    "\n",
+    "depths_uzet = uzdat_uzet.columns\n",
+    "depths_uzet = [\n",
+    "    round(float(d.strip().split(\"=\")[-1]), 2) for d in depths_uzet if \"=\" in d\n",
+    "]\n",
+    "\n",
+    "depths_theta = uzdat_theta.columns\n",
+    "depths_theta = [\n",
+    "    round(float(d.strip().split(\"=\")[-1]), 2) for d in depths_theta if \"=\" in d\n",
+    "]\n",
+    "\n",
+    "norm = plt.Normalize(0.5, 20.0)\n",
+    "\n",
+    "layelvs = [top] + botm\n",
+    "cell_half_thk = [\n",
+    "    round((up - down) / 2, 2) for (up, down) in zip(layelvs[:-1], layelvs[1:])\n",
+    "]\n",
+    "node_elvs = np.array(layelvs[:-1]) - np.array(cell_half_thk)\n",
+    "y = layelvs[0:20]\n",
+    "x = [0.5 for i in np.arange(len(y))]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "19bd33d3-09d4-4fea-a90d-b78cafa9acd8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Animate results\n",
+    "\n",
+    "matplotlib.rcParams[\"animation.embed_limit\"] = 2**128\n",
+    "\n",
+    "fig, (ax1, ax2, ax3, ax4) = plt.subplots(\n",
+    "    nrows=1,\n",
+    "    ncols=4,\n",
+    "    figsize=(8.5, 4),\n",
+    "    sharey=True,\n",
+    "    # gridspec_kw=dict(width_ratios=[3, 1, 1], wspace=0.2)\n",
+    ")\n",
+    "\n",
+    "# Look at the effects of Brooks-Corey epsilon\n",
+    "ax1.set_ylabel(\"Depth (m)\")\n",
+    "ax1.set_xlabel(\"Water Content\")\n",
+    "(l1,) = ax1.plot([], [], \"b-\", mfc=\"none\", label=r\"$\\theta; \\epsilon = 3.6$\")\n",
+    "(l2,) = ax1.plot([], [], \"r-\", mfc=\"none\", label=r\"$\\theta; \\epsilon = 7.2$\")\n",
+    "line_eps = [l1, l2]\n",
+    "ax1.legend(loc=\"lower right\")\n",
+    "ax1.set_xlim(0.0, 0.25)\n",
+    "ax1.set_ylim(52.0, 60.1)\n",
+    "for yline in layelvs[: len(y)]:\n",
+    "    ax1.axhline(yline, color=\"gray\", linewidth=0.5)\n",
+    "\n",
+    "\n",
+    "# Look at the effects of vertical hydraulic conductivity\n",
+    "ax2.set_xlabel(\"Water Content\")\n",
+    "(l3,) = ax2.plot(\n",
+    "    [], [], \"b-\", mfc=\"none\", label=r\"$\\theta; vks = 3e-6 \\dfrac{m}{s}$\"\n",
+    ")\n",
+    "(l4,) = ax2.plot(\n",
+    "    [], [], \"r-\", mfc=\"none\", label=r\"$\\theta; vks = 3e-5 \\dfrac{m}{s}$\"\n",
+    ")\n",
+    "line_vks = [l3, l4]\n",
+    "ax2.legend(loc=\"lower right\")\n",
+    "ax2.set_xlim(0.0, 0.25)\n",
+    "ax2.set_ylim(52.0, 60.1)\n",
+    "for yline in layelvs[: len(y)]:\n",
+    "    ax2.axhline(yline, color=\"gray\", linewidth=0.5)\n",
+    "\n",
+    "\n",
+    "# Look at the effects of ET from the unsaturated zone\n",
+    "ax3.set_xlabel(\"Water Content\")\n",
+    "(l5,) = ax3.plot([], [], \"b-\", mfc=\"none\", label=r\"$\\theta; No ET$\")\n",
+    "(l6,) = ax3.plot([], [], \"r-\", mfc=\"none\", label=r\"$\\theta; With ET$\")\n",
+    "line_uzet = [l5, l6]\n",
+    "ax3.legend(loc=\"lower right\")\n",
+    "ax3.set_xlim(0.0, 0.25)\n",
+    "ax3.set_ylim(52.0, 60.1)\n",
+    "for yline in layelvs[: len(y)]:\n",
+    "    ax3.axhline(yline, color=\"gray\", linewidth=0.5)\n",
+    "\n",
+    "\n",
+    "# Look at the effects of vertical hydraulic conductivity\n",
+    "ax4.set_xlabel(\"Water Content\")\n",
+    "(l7,) = ax4.plot([], [], \"b-\", mfc=\"none\", label=r\"$\\theta; \\theta_s = 0.25$\")\n",
+    "(l8,) = ax4.plot([], [], \"r-\", mfc=\"none\", label=r\"$\\theta; \\theta_s = 0.15$\")\n",
+    "line_theta = [l7, l8]\n",
+    "ax4.legend(loc=\"lower right\")\n",
+    "ax4.set_xlim(0.0, 0.25)\n",
+    "ax4.set_ylim(52.0, 60.1)\n",
+    "for yline in layelvs[: len(y)]:\n",
+    "    ax4.axhline(yline, color=\"gray\", linewidth=0.5)\n",
+    "\n",
+    "\n",
+    "times = uzdat_base.index\n",
+    "\n",
+    "\n",
+    "def init():\n",
+    "    ax1.set_title(\"Time = %.2f years\" % (times[0]))\n",
+    "\n",
+    "\n",
+    "def update(j):\n",
+    "    # Moisture Content\n",
+    "    XBase = uzdat_base.iloc[j].values\n",
+    "    XScen1 = uzdat_eps.iloc[j].values\n",
+    "    XScen2 = uzdat_vks.iloc[j].values\n",
+    "    XScen3 = uzdat_uzet.iloc[j].values\n",
+    "    XScen4 = uzdat_theta.iloc[j].values\n",
+    "    t = times[j] / 365\n",
+    "    line_eps[0].set_data(XBase, depths_base)\n",
+    "    line_eps[1].set_data(XScen1, depths_eps)\n",
+    "    ax1.set_title(\"Time = %.2f years\" % (t))\n",
+    "\n",
+    "    line_vks[0].set_data(XBase, depths_base)\n",
+    "    line_vks[1].set_data(XScen2, depths_vks)\n",
+    "    ax2.set_title(\"Time = %.2f years\" % (t))\n",
+    "\n",
+    "    line_uzet[0].set_data(XBase, depths_base)\n",
+    "    line_uzet[1].set_data(XScen3, depths_uzet)\n",
+    "    ax3.set_title(\"Time = %.2f years\" % (t))\n",
+    "\n",
+    "    line_theta[0].set_data(XBase, depths_base)\n",
+    "    line_theta[1].set_data(XScen4, depths_theta)\n",
+    "    ax4.set_title(\"Time = %.2f years\" % (t))\n",
+    "\n",
+    "    return\n",
+    "\n",
+    "\n",
+    "ani = animation.FuncAnimation(\n",
+    "    fig, update, interval=100, frames=365\n",
+    ")  # frames=len(times)\n",
+    "HTML(ani.to_jshtml())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "75d6ef7b-9595-464d-8c48-7af862cd4674",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ani.save(\"UZF6.gif\", writer='pillow',fps=60)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a64e8901-1574-4828-92ef-55453dee1a51",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/exercises/gwf_advanced/step1_uzf-mvr.ipynb b/exercises/gwf_advanced/step1_uzf-mvr.ipynb
new file mode 100644
index 0000000..aea4711
--- /dev/null
+++ b/exercises/gwf_advanced/step1_uzf-mvr.ipynb
@@ -0,0 +1,479 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "09cd4c0d-ed77-4b8d-9d65-fccdf103d008",
+   "metadata": {},
+   "source": [
+    "# U.S. Geological Survey Class GW3099\n",
+    "Advanced Modeling of Groundwater Flow (GW3099)\\\n",
+    "Boise, Idaho\\\n",
+    "September 16 - 20, 2024\n",
+    "\n",
+    "![title](../../images/ClassLocation.jpg)\n",
+    "\n",
+    "# Create a model with Advanced Packages - Part 2\n",
+    "\n",
+    "Load model from Part 1 and amend with UZF & MVR, rerun and look at the MVR budget"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ccf2fce9-eced-4882-a807-6b4719522412",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pathlib import Path\n",
+    "\n",
+    "import flopy\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3ee4cbe9-8b6e-4eb8-a807-ade0eae794a4",
+   "metadata": {},
+   "source": [
+    "## Load the model we just built"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "50d232a7-3398-45e0-ad59-3ee225a120cd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "name = \"ad-p1\"\n",
+    "ws = Path(f\"./temp/{name}\")\n",
+    "\n",
+    "sim = flopy.mf6.MFSimulation.load(\n",
+    "    sim_name=name,\n",
+    "    exe_name=\"mf6\",\n",
+    "    sim_ws=ws,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "39bb1ef2-26c4-4c9e-8855-c2058f15e671",
+   "metadata": {},
+   "source": [
+    "### Get the GWF model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f2e5cf5e-6aaa-49fc-9730-f3e836f5f813",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gwf = sim.get_model(\"ad-p1\")\n",
+    "\n",
+    "# get some needed parameter values\n",
+    "nlay = gwf.dis.nlay.array"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "64d87ffc-d8c1-4a52-ac35-835cb968fa11",
+   "metadata": {},
+   "source": [
+    "## Add instance of UZF package"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fa2928aa-0cbc-4221-8c18-8186838f4109",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Define which row/column index the stack of uzf objects will appear in.\n",
+    "# Recall that the pumping well is in row, column = 3, 13.\n",
+    "# Place an irrigated field one column to the right (downhill)\n",
+    "uzf_row, uzf_col = 3, 14\n",
+    "vks = 0.1  # 0.001\n",
+    "thts = gwf.sto.sy.array[0, 0, 0]\n",
+    "thtr = 0.01\n",
+    "thti = 0.05\n",
+    "eps = 4\n",
+    "\n",
+    "uzf_pkdat = []\n",
+    "\n",
+    "for iuzno in np.arange(0, nlay, 1):\n",
+    "    # only the top layer is needs land flag = 1\n",
+    "    if iuzno == 0:\n",
+    "        landflag = 1\n",
+    "        surfdep = 0.1\n",
+    "    else:\n",
+    "        landflag = 0\n",
+    "        surfdep = 0.0\n",
+    "\n",
+    "    # all but the bottom most layer have a vertical connection\n",
+    "    if iuzno < nlay - 1:\n",
+    "        ivertconn = iuzno + 1\n",
+    "    else:\n",
+    "        # bottom layer\n",
+    "        ivertconn = -1\n",
+    "\n",
+    "    # Add the UZF object to the packagedata list object\n",
+    "    uzf_pkdat.append(\n",
+    "        # ifno, cellid(l, r, c),           landflag,  ivertcon, surfdep, vks, thtr, thts, thti, eps\n",
+    "        [\n",
+    "            iuzno,\n",
+    "            (iuzno, uzf_row, uzf_col),\n",
+    "            landflag,\n",
+    "            ivertconn,\n",
+    "            surfdep,\n",
+    "            vks,\n",
+    "            thtr,\n",
+    "            thts,\n",
+    "            thti,\n",
+    "            eps,\n",
+    "        ]\n",
+    "    )\n",
+    "\n",
+    "# stress period data\n",
+    "finf = 0.0\n",
+    "pet = 0.0\n",
+    "extdp = 0.0\n",
+    "extwc = 0.0\n",
+    "zero = 0.0\n",
+    "\n",
+    "uzf_spd = {\n",
+    "    0: [[0, finf, pet, extdp, extwc, zero, zero, zero]],\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "50f24b0e-54a9-4101-bc20-d830222862bf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "uzf = flopy.mf6.ModflowGwfuzf(\n",
+    "    gwf,\n",
+    "    print_flows=True,\n",
+    "    save_flows=True,\n",
+    "    wc_filerecord=name + \".uzfwc.bin\",\n",
+    "    simulate_et=False,\n",
+    "    simulate_gwseep=False,\n",
+    "    linear_gwet=False,\n",
+    "    boundnames=False,\n",
+    "    ntrailwaves=15,\n",
+    "    nwavesets=40,\n",
+    "    # observations=uzf_obs,\n",
+    "    nuzfcells=len(uzf_pkdat),\n",
+    "    packagedata=uzf_pkdat,\n",
+    "    perioddata=uzf_spd,\n",
+    "    budget_filerecord=f\"{name}.uzf.bud\",\n",
+    "    pname=\"UZF\",\n",
+    "    filename=f\"{name}.uzf\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5b6e620c-6150-4bba-92dd-d52c140ccd56",
+   "metadata": {},
+   "source": [
+    "## Setup MVR package\n",
+    "\n",
+    "Retrieve package names:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "eca3319c-d25c-4479-a586-c2bebdab2121",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lak_nm = gwf.lak.package_name\n",
+    "sfr_nm = gwf.sfr.package_name\n",
+    "maw_nm = gwf.maw.package_name\n",
+    "uzf_nm = gwf.uzf.package_name"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "33a92159-c882-42e2-936c-958fb75962b6",
+   "metadata": {},
+   "source": [
+    "### Update packages included in the MVR package"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7181941b-4486-49ed-8102-8732cb636f75",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# activate the mover option within every package that will be used as a provider or a receiver\n",
+    "gwf.sfr.mover = True\n",
+    "gwf.lak.mover = True\n",
+    "gwf.maw.mover = True\n",
+    "gwf.uzf.mover = True\n",
+    "\n",
+    "# gwf.sfr.save_flows = True\n",
+    "\n",
+    "# add an outlet to the lake package\n",
+    "gwf.lak.noutlet = 1\n",
+    "#                  outletno, lakein, lakeout, couttype, invert, width, rough, slope\n",
+    "gwf.lak.outlets = [(0, 0, 0, \"manning\", 130.0, 50.0, 0.03, 0.01)]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3854d120-3672-470c-81af-58ec98b97753",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "packages = [(lak_nm,), (sfr_nm,), (maw_nm,), (uzf_nm,)]\n",
+    "pump_frac = 0.5\n",
+    "\n",
+    "mvr_pkdat = []\n",
+    "\n",
+    "# lak -> sfr (lake outlet)\n",
+    "mvr_pkdat.append([lak_nm, 0, sfr_nm, 0, \"factor\", 1.0])\n",
+    "\n",
+    "# maw -> uzf (gw irrigation)\n",
+    "mvr_pkdat.append([maw_nm, 0, uzf_nm, 0, \"factor\", pump_frac])\n",
+    "\n",
+    "# maw -> sfr (pump to stream for delivery somewhere downstream)\n",
+    "mvr_pkdat.append([maw_nm, 0, sfr_nm, 2, \"factor\", 1 - pump_frac])\n",
+    "\n",
+    "# uzf -> sfr (irrigation tail-water plus any other runoff\n",
+    "mvr_pkdat.append([uzf_nm, 0, sfr_nm, 2, \"factor\", 1.0])\n",
+    "\n",
+    "flopy.mf6.ModflowGwfmvr(\n",
+    "    gwf,\n",
+    "    maxmvr=len(mvr_pkdat),\n",
+    "    budget_filerecord=f\"{name}.mvr.bud\",\n",
+    "    maxpackages=len(packages),\n",
+    "    print_flows=True,\n",
+    "    packages=packages,\n",
+    "    perioddata=mvr_pkdat,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3a871bef-ab79-4caf-b7af-5c33a7736419",
+   "metadata": {},
+   "source": [
+    "## Write the simulation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fd2eb7f0-ff14-44e7-b706-8d12a083a485",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim.write_simulation()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e8dc3726-3a38-4a7a-995a-750df45d7d7b",
+   "metadata": {},
+   "source": [
+    "## Run the simulation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fc541f49-76c5-4c89-98cf-ed4a7e43c51c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim.run_simulation()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "954585bd-c65d-4336-87ec-42ffbeffad39",
+   "metadata": {},
+   "source": [
+    "## Have another look at the output"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "20c402f3-6606-416f-bc18-6c55bc78dece",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "extents = (0.0, gwf.dis.delr.array.sum(), 0.0, gwf.dis.delc.array.sum())\n",
+    "\n",
+    "# load the observations\n",
+    "lak_results = gwf.lak.output.obs().data\n",
+    "sfr_results = gwf.sfr.output.obs().data\n",
+    "gwf_results = gwf.obs[0].output.obs().data\n",
+    "\n",
+    "# Figure properties\n",
+    "figure_size = (6.3, 5.6)\n",
+    "masked_values = (0, 1e30, -1e30)\n",
+    "\n",
+    "# create MODFLOW 6 head object\n",
+    "hobj = gwf.output.head()\n",
+    "\n",
+    "# create MODFLOW 6 cell-by-cell budget object\n",
+    "cobj = gwf.output.budget()\n",
+    "\n",
+    "kstpkper = hobj.get_kstpkper()\n",
+    "\n",
+    "head = hobj.get_data(kstpkper=kstpkper[0])\n",
+    "qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(\n",
+    "    cobj.get_data(text=\"DATA-SPDIS\", kstpkper=kstpkper[0])[0],\n",
+    "    gwf,\n",
+    ")\n",
+    "\n",
+    "# add lake stage to heads\n",
+    "head[head == 1e30] = lak_results[\"STAGE\"][-1]\n",
+    "\n",
+    "# observation locations\n",
+    "xcenters, ycenters = gwf.modelgrid.xycenters[0], gwf.modelgrid.xycenters[1]\n",
+    "p1 = (xcenters[3], ycenters[3])\n",
+    "p2 = (xcenters[13], ycenters[13])\n",
+    "\n",
+    "shape3d = (gwf.dis.nlay.array, gwf.dis.nrow.array, gwf.dis.ncol.array)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2a46f3fb-9cc9-456c-8d60-61bfbb3b83da",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axd = plt.subplot_mosaic(\n",
+    "    [\n",
+    "        [\"a\"],\n",
+    "        [\"a\"],\n",
+    "        [\"b\"],\n",
+    "    ],\n",
+    "    layout=\"constrained\",\n",
+    "    figsize=(4, 6.9),\n",
+    ")\n",
+    "\n",
+    "ax = axd[\"a\"]\n",
+    "mm = flopy.plot.PlotMapView(gwf, ax=ax, extent=extents)\n",
+    "mm.plot_bc(\"CHD\", color=\"cyan\")\n",
+    "mm.plot_bc(\"SFR\", color=\"blue\", alpha=0.1)\n",
+    "mm.plot_bc(\"UZF\", color=\"lightgreen\", alpha=0.7)\n",
+    "mm.plot_inactive(color_noflow=\"#5DBB63\")\n",
+    "mm.plot_grid(lw=0.5, color=\"black\")\n",
+    "cv = mm.contour_array(\n",
+    "    head,\n",
+    "    levels=np.arange(140, 160, 2),\n",
+    "    linewidths=0.75,\n",
+    "    linestyles=\"-\",\n",
+    "    colors=\"blue\",\n",
+    ")\n",
+    "plt.clabel(cv, fmt=\"%1.0f\")\n",
+    "mm.plot_vector(qx, qy, normalize=True, color=\"0.75\")\n",
+    "ax.plot(p1[0], p1[1], marker=\"o\", mfc=\"red\", mec=\"black\", ms=4)\n",
+    "ax.annotate(\"Point A\", (p1[0] + 150, p1[1]))\n",
+    "ax.plot(p2[0], p2[1], marker=\"o\", mfc=\"red\", mec=\"black\", ms=4)\n",
+    "ax.annotate(\"Point B\", (p2[0] + 150, p2[1]))\n",
+    "ax.plot(p2[0], p1[1], marker=\"o\", mfc=\"yellow\", mec=\"purple\", ms=4)\n",
+    "ax.annotate(\"MAW\", (p2[0] - 1500, p1[1] + 150))\n",
+    "ax.annotate(\"UZF\", (10.75e3, 11e3))\n",
+    "ax.set_xlabel(\"x-coordinate, in feet\")\n",
+    "ax.set_ylabel(\"y-coordinate, in feet\")\n",
+    "\n",
+    "ax = axd[\"b\"]\n",
+    "xs = flopy.plot.PlotCrossSection(gwf, ax=ax, line={\"row\": 8})\n",
+    "xs.plot_array(np.ones(shape3d), head=head, cmap=\"jet\")\n",
+    "xs.plot_bc(\"CHD\", color=\"cyan\", head=head)\n",
+    "xs.plot_ibound(color_noflow=\"#5DBB63\", head=head)\n",
+    "xs.plot_grid(lw=0.5, color=\"black\")\n",
+    "ax.set_xlabel(\"x-coordinate, in feet\")\n",
+    "ax.set_ylim(67, 160)\n",
+    "ax.set_ylabel(\"Elevation, in feet\")\n",
+    "\n",
+    "plt.show(block=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d66867d2-428a-428e-b2a6-a4b14d64de69",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(sim.sim_path)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c22dc53c-ad5e-4dc9-8eec-3a1762b4c37a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lst_pth = ws / f\"{name}.lst\"\n",
+    "budget_key = \"WATER MOVER BUDGET FOR ENTIRE MODEL AT END OF TIME STEP  100\"\n",
+    "lst = flopy.utils.MfListBudget(lst_pth, budgetkey=budget_key)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c5446a1a-cea0-4067-be36-3aa98501a78a",
+   "metadata": {},
+   "source": [
+    "### Budget for current stress period"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9501e034-adc1-4c5f-b6f3-a64195f5b5f7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lst.get_dataframes()[0]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c2b47b6b-e3cd-44a7-a2bd-5df2a6a37c3d",
+   "metadata": {},
+   "source": [
+    "### Cumulative budget to the current stress period\n",
+    "(which happens to be the last stress period)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "efb3c579-20cb-4a62-984f-39c37a07e3f3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lst.get_dataframes()[1]"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/exercises/gwf_advanced/step4_more-sfr.ipynb b/exercises/gwf_advanced/step4_more-sfr.ipynb
new file mode 100644
index 0000000..9bc030e
--- /dev/null
+++ b/exercises/gwf_advanced/step4_more-sfr.ipynb
@@ -0,0 +1,542 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "09cd4c0d-ed77-4b8d-9d65-fccdf103d008",
+   "metadata": {},
+   "source": [
+    "# U.S. Geological Survey Class GW3099\n",
+    "Advanced Modeling of Groundwater Flow (GW3099)\\\n",
+    "Boise, Idaho\\\n",
+    "September 16 - 20, 2024\n",
+    "\n",
+    "![title](../../images/ClassLocation.jpg)\n",
+    "\n",
+    "# Create a model with Advanced Packages - Part 2\n",
+    "\n",
+    "Load model from Part 1 and amend with UZF & MVR, rerun and look at the MVR budget"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ccf2fce9-eced-4882-a807-6b4719522412",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import warnings\n",
+    "from pathlib import Path\n",
+    "\n",
+    "import flopy\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "warnings.simplefilter(\"ignore\", UserWarning)\n",
+    "warnings.simplefilter(\"ignore\", DeprecationWarning)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3ee4cbe9-8b6e-4eb8-a807-ade0eae794a4",
+   "metadata": {},
+   "source": [
+    "## Load the model amended with UZF and MVR"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "50d232a7-3398-45e0-ad59-3ee225a120cd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "name = \"ad-p1\"\n",
+    "ws = Path(f\"./temp/{name}\")\n",
+    "\n",
+    "sim = flopy.mf6.MFSimulation.load(\n",
+    "    sim_name=name,\n",
+    "    exe_name=\"mf6\",\n",
+    "    sim_ws=ws,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "39bb1ef2-26c4-4c9e-8855-c2058f15e671",
+   "metadata": {},
+   "source": [
+    "### Get the GWF model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f2e5cf5e-6aaa-49fc-9730-f3e836f5f813",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gwf = sim.get_model(\"ad-p1\")\n",
+    "\n",
+    "# get some needed parameter values\n",
+    "nlay = gwf.dis.nlay.array"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ff341bd6-9c91-4462-8b83-d7c8f016d0c8",
+   "metadata": {},
+   "source": [
+    "### Add a diversion within the SFR package\n",
+    "For simplicity, diversion will run due north away from the existing stream that drains the lake.  The diversion will traverse 4 cells all in column 14 (0-based)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d9b291bc-03f3-4055-a326-dab359e8c46a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Original data\n",
+    "bed_elev = [142.0, 141.5, 141.0, 140.5, 140.0] + [140.4, 140.3, 140.2, 140.1]\n",
+    "upstfr = [0.0, 1.0, 1.0, 1.0, 1.0] + [0.0, 1.0, 1.0, 1.0]\n",
+    "nconn = [1, 2, 2, 3, 1] + [2, 2, 2, 1]\n",
+    "ndiv = [0, 0, 0, 1, 0] + [0, 0, 0, 0]\n",
+    "nreaches = len(bed_elev)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "adeea311-e615-49ea-9acd-94e898348b60",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# <ifno> <cellid(ncelldim)> <rlen> <rwid> <rgrd> <rtp> <rbth> <rhk> <man> <ncon> <ustrf> <ndv>\n",
+    "sfr_pakdata = []\n",
+    "ct = 0\n",
+    "divCol = 14  # 0-based\n",
+    "for idx in range(nreaches):\n",
+    "    rchlen = 1000.0\n",
+    "    if idx <= 4:\n",
+    "        cellid = (0, 8, 11 + idx)\n",
+    "    else:\n",
+    "        ct += 1\n",
+    "        cellid = (0, 8 - ct, divCol)\n",
+    "        if idx <= 6:\n",
+    "            rchlen = 500.0\n",
+    "\n",
+    "    sfr_pakdata.append(\n",
+    "        (\n",
+    "            idx,\n",
+    "            cellid,\n",
+    "            rchlen,\n",
+    "            1.0,\n",
+    "            1e-3,\n",
+    "            bed_elev[idx],\n",
+    "            0.1,\n",
+    "            1.0,\n",
+    "            0.035,\n",
+    "            nconn[idx],\n",
+    "            upstfr[idx],\n",
+    "            ndiv[idx],\n",
+    "            \"sfr\",\n",
+    "        )\n",
+    "    )\n",
+    "\n",
+    "sfr_pakdata"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bfc38ce6-c83f-4358-b76a-a7cb121abf30",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sfr_conn = [\n",
+    "    [0, -1],\n",
+    "    [1, 0, -2],\n",
+    "    [2, 1, -3],\n",
+    "    [3, 2, -4, -5],\n",
+    "    [4, 3],\n",
+    "    [5, 3, -6],\n",
+    "    [6, 5, -7],\n",
+    "    [7, 6, -8],\n",
+    "    [8, 7],\n",
+    "]\n",
+    "sfr_conn"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "efd53ce0-d11d-4f09-b2b0-2ff3272104a6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "obs_file = f\"{name}.sfr.obs\"\n",
+    "csv_file = obs_file + \".csv\"\n",
+    "obs_dict = {\n",
+    "    obs_file + \".csv\": [\n",
+    "        (\"inflow_rch4\", \"inflow\", (3,)),\n",
+    "        (\"outflow_rch4\", \"outflow\", (3,)),\n",
+    "        (\"inflow_rch5\", \"inflow\", (4,)),\n",
+    "        (\"divertedAmt\", \"inflow\", (5,)),\n",
+    "    ]\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "02164c8f-a009-468d-81f3-02f14d76e0bb",
+   "metadata": {},
+   "source": [
+    "#### Define the `DIVERSIONS` block\n",
+    "Use the `CPRIOR` option `EXCESS`: <br>\n",
+    "<br>\n",
+    "`EXCESS`: a diversion is made only if $Q_{DS}$ for reach `IFNO` exceeds the value of `DIVFLOW`. If this occurs, then the quantity of water diverted is\n",
+    "the excess flow ($Q_{DS}$ − DIVFLOW) and $Q_{DS}$ from reach IFNO is set equal to DIVFLOW. This represents a flood-control type of diversion"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ba451d66-ec0c-482c-a756-e8fec35c5b8c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#      <ifno> <idv> <iconr> <cprior>\n",
+    "divs = [3, 0, 5, \"EXCESS\"]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9008e022-6428-4aae-82f7-a5c6efd16b19",
+   "metadata": {},
+   "source": [
+    "#### Define the stress period data\n",
+    "Because we added a diversion, need to set the transient data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c2e23628-681a-4e05-b6b1-a132210d87c9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "excess_amt = 10000.0\n",
+    "sfr_spd = {0: [3, \"DIVERSION\", 0, excess_amt]}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1b66ef3d-2ef4-42df-91ec-674c36a91fdc",
+   "metadata": {},
+   "source": [
+    "#### Redefine the SFR package\n",
+    "There is one gotcha here.  In the script that follows, if we don't specify the argument `pname`, FloPy will append the package to the existing `MODFLOW 6` simulation rather than replace the existing SFR package with the new one that includes the diversion."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "44504805-1d9c-41e2-91f3-004972c08b51",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sfr = flopy.mf6.ModflowGwfsfr(\n",
+    "    gwf,\n",
+    "    mover=True,\n",
+    "    boundnames=True,\n",
+    "    print_input=True,\n",
+    "    print_flows=True,\n",
+    "    print_stage=True,\n",
+    "    length_conversion=3.28081,\n",
+    "    time_conversion=86400.0,\n",
+    "    nreaches=nreaches,\n",
+    "    packagedata=sfr_pakdata,\n",
+    "    connectiondata=sfr_conn,\n",
+    "    diversions=divs,\n",
+    "    observations=obs_dict,\n",
+    "    perioddata=sfr_spd,\n",
+    "    pname=\"sfr_0\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3a871bef-ab79-4caf-b7af-5c33a7736419",
+   "metadata": {},
+   "source": [
+    "## Write the simulation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fd2eb7f0-ff14-44e7-b706-8d12a083a485",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim.write_simulation()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e8dc3726-3a38-4a7a-995a-750df45d7d7b",
+   "metadata": {},
+   "source": [
+    "## Run the simulation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fc541f49-76c5-4c89-98cf-ed4a7e43c51c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim.run_simulation()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "954585bd-c65d-4336-87ec-42ffbeffad39",
+   "metadata": {},
+   "source": [
+    "## Have another look at the output"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "20c402f3-6606-416f-bc18-6c55bc78dece",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "extents = (0.0, gwf.dis.delr.array.sum(), 0.0, gwf.dis.delc.array.sum())\n",
+    "\n",
+    "# load the observations\n",
+    "lak_results = gwf.lak.output.obs().data\n",
+    "sfr_results = gwf.sfr.output.obs().data\n",
+    "gwf_results = gwf.obs[0].output.obs().data\n",
+    "\n",
+    "# Figure properties\n",
+    "figure_size = (6.3, 5.6)\n",
+    "masked_values = (0, 1e30, -1e30)\n",
+    "\n",
+    "# create MODFLOW 6 head object\n",
+    "hobj = gwf.output.head()\n",
+    "\n",
+    "# create MODFLOW 6 cell-by-cell budget object\n",
+    "cobj = gwf.output.budget()\n",
+    "\n",
+    "kstpkper = hobj.get_kstpkper()\n",
+    "\n",
+    "head = hobj.get_data(kstpkper=kstpkper[0])\n",
+    "qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(\n",
+    "    cobj.get_data(text=\"DATA-SPDIS\", kstpkper=kstpkper[0])[0],\n",
+    "    gwf,\n",
+    ")\n",
+    "\n",
+    "# add lake stage to heads\n",
+    "head[head == 1e30] = lak_results[\"STAGE\"][-1]\n",
+    "\n",
+    "# observation locations\n",
+    "xcenters, ycenters = gwf.modelgrid.xycenters[0], gwf.modelgrid.xycenters[1]\n",
+    "p1 = (xcenters[3], ycenters[3])\n",
+    "p2 = (xcenters[13], ycenters[13])\n",
+    "\n",
+    "shape3d = (gwf.dis.nlay.array, gwf.dis.nrow.array, gwf.dis.ncol.array)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2a46f3fb-9cc9-456c-8d60-61bfbb3b83da",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axd = plt.subplot_mosaic(\n",
+    "    [\n",
+    "        [\"a\"],\n",
+    "        [\"a\"],\n",
+    "        [\"b\"],\n",
+    "    ],\n",
+    "    layout=\"constrained\",\n",
+    "    figsize=(4, 6.9),\n",
+    ")\n",
+    "\n",
+    "ax = axd[\"a\"]\n",
+    "mm = flopy.plot.PlotMapView(gwf, ax=ax, extent=extents)\n",
+    "mm.plot_bc(\"CHD\", color=\"cyan\")\n",
+    "mm.plot_bc(\"SFR\", color=\"blue\", alpha=0.1)\n",
+    "mm.plot_bc(\"UZF\", color=\"lightgreen\", alpha=0.7)\n",
+    "mm.plot_inactive(color_noflow=\"#5DBB63\")\n",
+    "mm.plot_grid(lw=0.5, color=\"black\")\n",
+    "cv = mm.contour_array(\n",
+    "    head,\n",
+    "    levels=np.arange(140, 160, 2),\n",
+    "    linewidths=0.75,\n",
+    "    linestyles=\"-\",\n",
+    "    colors=\"blue\",\n",
+    ")\n",
+    "plt.clabel(cv, fmt=\"%1.0f\")\n",
+    "mm.plot_vector(qx, qy, normalize=True, color=\"0.75\")\n",
+    "ax.plot(p1[0], p1[1], marker=\"o\", mfc=\"red\", mec=\"black\", ms=4)\n",
+    "ax.annotate(\"Point A\", (p1[0] + 150, p1[1]))\n",
+    "ax.plot(p2[0], p2[1], marker=\"o\", mfc=\"red\", mec=\"black\", ms=4)\n",
+    "ax.annotate(\"Point B\", (p2[0] + 150, p2[1]))\n",
+    "ax.plot(p2[0], p1[1], marker=\"o\", mfc=\"yellow\", mec=\"purple\", ms=4)\n",
+    "ax.annotate(\"MAW\", (p2[0] - 1500, p1[1] + 150))\n",
+    "ax.annotate(\"UZF\", (10.75e3, 11e3))\n",
+    "ax.set_xlabel(\"x-coordinate, in feet\")\n",
+    "ax.set_ylabel(\"y-coordinate, in feet\")\n",
+    "\n",
+    "ax = axd[\"b\"]\n",
+    "xs = flopy.plot.PlotCrossSection(gwf, ax=ax, line={\"row\": 8})\n",
+    "xs.plot_array(np.ones(shape3d), head=head, cmap=\"jet\")\n",
+    "xs.plot_bc(\"CHD\", color=\"cyan\", head=head)\n",
+    "xs.plot_ibound(color_noflow=\"#5DBB63\", head=head)\n",
+    "xs.plot_grid(lw=0.5, color=\"black\")\n",
+    "ax.set_xlabel(\"x-coordinate, in feet\")\n",
+    "ax.set_ylim(67, 160)\n",
+    "ax.set_ylabel(\"Elevation, in feet\")\n",
+    "\n",
+    "plt.show(block=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f0b3c713-31bd-4c59-80ce-672d8685a981",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sfr_results = gwf.sfr.output.obs().data\n",
+    "sfr_results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "886d7a91-9331-4dd4-8f08-7ff447e63575",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dtype = [\n",
+    "    (\"time\", float),\n",
+    "    (\"INFLOW_RCH4\", float),\n",
+    "    (\"OUTFLOW_RCH4\", float),\n",
+    "    (\"INFLOW_RCH5\", float),\n",
+    "    (\"DIVERTEDAMT\", float),\n",
+    "]\n",
+    "\n",
+    "results = np.zeros((sfr_results.shape[0] + 1), dtype=dtype)\n",
+    "results[\"time\"][1:] = sfr_results[\"totim\"]\n",
+    "results[\"INFLOW_RCH4\"][1:] = sfr_results[\"INFLOW_RCH4\"]\n",
+    "results[\"OUTFLOW_RCH4\"][1:] = sfr_results[\"OUTFLOW_RCH4\"]\n",
+    "results[\"INFLOW_RCH5\"][1:] = sfr_results[\"INFLOW_RCH5\"]\n",
+    "results[\"DIVERTEDAMT\"][1:] = sfr_results[\"DIVERTEDAMT\"]\n",
+    "results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a40e5356-51bc-4c9c-acf3-78aa152e87b1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create the figure\n",
+    "fig, ax = plt.subplots(\n",
+    "    ncols=1,\n",
+    "    nrows=1,\n",
+    "    sharex=True,\n",
+    "    figsize=(6.3, 3.15),\n",
+    "    constrained_layout=True,\n",
+    ")\n",
+    "\n",
+    "ax.set_xlim(0, 3000)\n",
+    "# ax.set_ylim(110, 160)\n",
+    "ax.plot(\n",
+    "    results[\"time\"],\n",
+    "    results[\"INFLOW_RCH4\"],\n",
+    "    lw=0.75,\n",
+    "    ls=\"--\",\n",
+    "    color=\"black\",\n",
+    "    label=\"Reach 4 inflow\",\n",
+    ")\n",
+    "ax.plot(\n",
+    "    results[\"time\"],\n",
+    "    results[\"OUTFLOW_RCH4\"],\n",
+    "    lw=0.75,\n",
+    "    ls=\":\",\n",
+    "    color=\"red\",\n",
+    "    label=\"Reach 4 outflow\",\n",
+    ")\n",
+    "ax.plot(\n",
+    "    results[\"time\"],\n",
+    "    results[\"INFLOW_RCH5\"],\n",
+    "    lw=1.0,\n",
+    "    ls=\"-\",\n",
+    "    color=\"green\",\n",
+    "    label=\"Reach 5 inflow\",\n",
+    ")\n",
+    "ax.plot(\n",
+    "    results[\"time\"],\n",
+    "    results[\"DIVERTEDAMT\"],\n",
+    "    lw=1.0,\n",
+    "    ls=\"-\",\n",
+    "    color=\"blue\",\n",
+    "    label=\"Diversion\",\n",
+    ")\n",
+    "ax.axhline(y=0, color=\"black\")\n",
+    "ax.set_xlabel(\"Simulation time, in days\")\n",
+    "ax.set_ylabel(\"Flow, in $ft^3$ per day\")\n",
+    "plt.legend(loc=\"lower left\")\n",
+    "\n",
+    "plt.show(block=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "33a92159-c882-42e2-936c-958fb75962b6",
+   "metadata": {},
+   "source": [
+    "### Update MVR package"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b8a853a6-161c-4838-adda-d7fddec88e38",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}