From c66e60fc5c098cf9ef1ee605120d7fd0bf9ef197 Mon Sep 17 00:00:00 2001
From: Daniel Buckley <buckley.daniel.n@gmail.com>
Date: Sat, 13 Jan 2024 16:17:04 -0700
Subject: [PATCH] Implement Newton-Raphson version of BCF2SS

Fixed OC file extension
---
 .../gwf_bcf2ss/gwf_bcf2ss_p1a.ipynb           |  17 +-
 .../gwf_bcf2ss/gwf_bcf2ss_p2a.ipynb           | 805 ++++++++++++++++++
 .../gwf_bcf2ss/model_p1a/gwf_bcf2ss.dis       |   2 +-
 .../gwf_bcf2ss/model_p1a/gwf_bcf2ss.ic        |   2 +-
 .../gwf_bcf2ss/model_p1a/gwf_bcf2ss.ims       |   2 +-
 .../gwf_bcf2ss/model_p1a/gwf_bcf2ss.nam       |  14 +-
 .../gwf_bcf2ss/model_p1a/gwf_bcf2ss.npf       |   2 +-
 .../model_p1a/{gwf_bcf2ss => gwf_bcf2ss.oc}   |   2 +-
 .../gwf_bcf2ss/model_p1a/gwf_bcf2ss.rcha      |   2 +-
 .../gwf_bcf2ss/model_p1a/gwf_bcf2ss.riv       |   2 +-
 .../gwf_bcf2ss/model_p1a/gwf_bcf2ss.tdis      |   2 +-
 .../gwf_bcf2ss/model_p1a/gwf_bcf2ss.wel       |   2 +-
 .../gwf_bcf2ss/model_p1a/mfsim.nam            |   2 +-
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.dis       |  23 +
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.dis.grb   | Bin 0 -> 16048 bytes
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.ic        |   6 +
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.ims       |  17 +
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.lst       | 377 ++++++++
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.nam       |  16 +
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.npf       |  16 +
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.oc        |  11 +
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.rcha      |  10 +
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.riv       |  21 +
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.tdis      |  14 +
 .../gwf_bcf2ss/model_p2a/gwf_bcf2ss.wel       |  13 +
 .../gwf_bcf2ss/model_p2a/mfsim.lst            | 234 +++++
 .../gwf_bcf2ss/model_p2a/mfsim.nam            |  19 +
 notebooks/learning/gwf_bcf2ss.ipynb           | 168 +++-
 28 files changed, 1778 insertions(+), 23 deletions(-)
 create mode 100644 notebooks/implementing/gwf_bcf2ss/gwf_bcf2ss_p2a.ipynb
 rename notebooks/implementing/gwf_bcf2ss/model_p1a/{gwf_bcf2ss => gwf_bcf2ss.oc} (70%)
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.dis
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.dis.grb
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.ic
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.ims
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.lst
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.nam
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.npf
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.oc
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.rcha
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.riv
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.tdis
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.wel
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/mfsim.lst
 create mode 100644 notebooks/implementing/gwf_bcf2ss/model_p2a/mfsim.nam

diff --git a/notebooks/implementing/gwf_bcf2ss/gwf_bcf2ss_p1a.ipynb b/notebooks/implementing/gwf_bcf2ss/gwf_bcf2ss_p1a.ipynb
index 6fc7203..32aecc6 100644
--- a/notebooks/implementing/gwf_bcf2ss/gwf_bcf2ss_p1a.ipynb
+++ b/notebooks/implementing/gwf_bcf2ss/gwf_bcf2ss_p1a.ipynb
@@ -625,14 +625,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "WARNING: Package with type oc already exists. Replacing existing package.\n"
+     ]
+    },
     {
      "data": {
       "text/plain": [
        "package_name = oc\n",
-       "filename = gwf_bcf2ss\n",
+       "filename = gwf_bcf2ss.oc\n",
        "package_type = oc\n",
        "model_or_simulation_package = model\n",
        "model_name = gwf_bcf2ss\n",
@@ -661,7 +668,7 @@
        "None\n"
       ]
      },
-     "execution_count": 15,
+     "execution_count": 22,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -675,7 +682,7 @@
     "        (\"head\", \"all\"),\n",
     "        (\"budget\", \"all\")\n",
     "    ],\n",
-    "    filename=f\"{sim_name}\",\n",
+    "    filename=f\"{sim_name}.oc\",\n",
     "    pname=\"oc\"\n",
     ")"
    ]
@@ -727,7 +734,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 23,
    "metadata": {},
    "outputs": [
     {
diff --git a/notebooks/implementing/gwf_bcf2ss/gwf_bcf2ss_p2a.ipynb b/notebooks/implementing/gwf_bcf2ss/gwf_bcf2ss_p2a.ipynb
new file mode 100644
index 0000000..8713bc9
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/gwf_bcf2ss_p2a.ipynb
@@ -0,0 +1,805 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from os import getcwd\n",
+    "from pathlib import Path\n",
+    "\n",
+    "import flopy as fp"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "WindowsPath('c:/Users/buckl/Documents/Github/flopy_mf6_work/notebooks/implementing/gwf_bcf2ss/model_p2a')"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "workspace = Path(getcwd(), \"model_p2a\")\n",
+    "workspace"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "sim_name = gwf_bcf2ss\n",
+       "sim_path = c:\\Users\\buckl\\Documents\\Github\\flopy_mf6_work\\notebooks\\implementing\\gwf_bcf2ss\\model_p2a\n",
+       "exe_name = C:\\Users\\buckl\\Documents\\Software\\mf6.4.1\\bin\\mf6.exe\n",
+       "\n",
+       "###################\n",
+       "Package mfsim.nam\n",
+       "###################\n",
+       "\n",
+       "package_name = mfsim.nam\n",
+       "filename = mfsim.nam\n",
+       "package_type = nam\n",
+       "model_or_simulation_package = simulation\n",
+       "simulation_name = gwf_bcf2ss\n",
+       "\n"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "sim_name = \"gwf_bcf2ss\"\n",
+    "sim = fp.mf6.MFSimulation(\n",
+    "    sim_name=sim_name,\n",
+    "    version=\"mf6\",\n",
+    "    exe_name=Path(\"C:\\\\Users\\\\buckl\\\\Documents\\\\Software\\\\mf6.4.1\\\\bin\\\\mf6.exe\"),\n",
+    "    sim_ws=workspace\n",
+    ")\n",
+    "sim"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "package_name = tdis\n",
+       "filename = gwf_bcf2ss.tdis\n",
+       "package_type = tdis\n",
+       "model_or_simulation_package = simulation\n",
+       "simulation_name = gwf_bcf2ss\n",
+       "\n",
+       "Block options\n",
+       "--------------------\n",
+       "time_units\n",
+       "{internal}\n",
+       "('days')\n",
+       "\n",
+       "\n",
+       "Block dimensions\n",
+       "--------------------\n",
+       "nper\n",
+       "{internal}\n",
+       "(2)\n",
+       "\n",
+       "\n",
+       "Block perioddata\n",
+       "--------------------\n",
+       "perioddata\n",
+       "{internal}\n",
+       "(rec.array([(1., 1, 1.), (1., 1, 1.)],\n",
+       "          dtype=[('perlen', '<f8'), ('nstp', '<i4'), ('tsmult', '<f8')]))\n",
+       "\n"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fp.mf6.ModflowTdis(\n",
+    "    simulation=sim,\n",
+    "    time_units=\"days\",\n",
+    "    nper=2,\n",
+    "    perioddata=[\n",
+    "        (1.0, 1, 1.0),\n",
+    "        (1.0, 1, 1.0)\n",
+    "    ],\n",
+    "    filename=f\"{sim_name}.tdis\",\n",
+    "    pname=\"tdis\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "package_name = ims\n",
+       "filename = gwf_bcf2ss.ims\n",
+       "package_type = ims\n",
+       "model_or_simulation_package = simulation\n",
+       "simulation_name = gwf_bcf2ss\n",
+       "\n",
+       "Block nonlinear\n",
+       "--------------------\n",
+       "outer_dvclose\n",
+       "{internal}\n",
+       "(1e-06)\n",
+       "\n",
+       "outer_maximum\n",
+       "{internal}\n",
+       "(500)\n",
+       "\n",
+       "\n",
+       "Block linear\n",
+       "--------------------\n",
+       "inner_maximum\n",
+       "{internal}\n",
+       "(100)\n",
+       "\n",
+       "inner_dvclose\n",
+       "{internal}\n",
+       "(1e-06)\n",
+       "\n",
+       "rcloserecord\n",
+       "{internal}\n",
+       "(rec.array([('inner_rclose', 0.001, 'strict')],\n",
+       "          dtype=[('inner_rclose_label', 'O'), ('inner_rclose', '<f8'), ('rclose_option', 'O')]))\n",
+       "\n",
+       "linear_acceleration\n",
+       "{internal}\n",
+       "('bicgstab')\n",
+       "\n",
+       "relaxation_factor\n",
+       "{internal}\n",
+       "(0.97)\n",
+       "\n"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fp.mf6.ModflowIms(\n",
+    "    simulation=sim,\n",
+    "    outer_dvclose=1e-6,\n",
+    "    outer_maximum=500,\n",
+    "    inner_maximum=100,\n",
+    "    inner_dvclose=1e-6,\n",
+    "    rcloserecord=[0.001, \"strict\"],\n",
+    "    linear_acceleration=\"bicgstab\",\n",
+    "    relaxation_factor=0.97,\n",
+    "    filename=f\"{sim_name}.ims\",\n",
+    "    pname=\"ims\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "name = gwf_bcf2ss\n",
+       "model_type = gwf6\n",
+       "version = mf6\n",
+       "model_relative_path = .\n"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "model = fp.mf6.ModflowGwf(\n",
+    "    simulation=sim,\n",
+    "    modelname=sim_name,\n",
+    "    model_nam_file=f\"{sim_name}.nam\",\n",
+    "    # from original nam file\n",
+    "    save_flows=True,\n",
+    "    # from original nam file\n",
+    "    newtonoptions=\"NEWTON\"\n",
+    ")\n",
+    "model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "package_name = dis\n",
+       "filename = gwf_bcf2ss.dis\n",
+       "package_type = dis\n",
+       "model_or_simulation_package = model\n",
+       "model_name = gwf_bcf2ss\n",
+       "\n",
+       "Block options\n",
+       "--------------------\n",
+       "length_units\n",
+       "{internal}\n",
+       "('feet')\n",
+       "\n",
+       "\n",
+       "Block dimensions\n",
+       "--------------------\n",
+       "nlay\n",
+       "{internal}\n",
+       "(2)\n",
+       "\n",
+       "nrow\n",
+       "{internal}\n",
+       "(10)\n",
+       "\n",
+       "ncol\n",
+       "{internal}\n",
+       "(15)\n",
+       "\n",
+       "\n",
+       "Block griddata\n",
+       "--------------------\n",
+       "delr\n",
+       "{constant 500.0}\n",
+       "\n",
+       "delc\n",
+       "{constant 500.0}\n",
+       "\n",
+       "top\n",
+       "{constant 150.0}\n",
+       "\n",
+       "botm\n",
+       "Layer_1{constant 50.0}\n",
+       "Layer_2{constant -50.0}\n",
+       "\n"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fp.mf6.ModflowGwfdis(\n",
+    "    model=model,\n",
+    "    length_units=\"feet\",\n",
+    "    nlay=2,\n",
+    "    nrow=10,\n",
+    "    ncol=15,\n",
+    "    delr=500.0,\n",
+    "    delc=500.0,\n",
+    "    top=150.0,\n",
+    "    botm = [\n",
+    "        50.0,\n",
+    "        -50.0\n",
+    "    ],\n",
+    "    filename=f\"{sim_name}.dis\",\n",
+    "    pname=\"dis\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "package_name = npf\n",
+       "filename = gwf_bcf2ss.npf\n",
+       "package_type = npf\n",
+       "model_or_simulation_package = model\n",
+       "model_name = gwf_bcf2ss\n",
+       "\n",
+       "Block options\n",
+       "--------------------\n",
+       "save_specific_discharge\n",
+       "{internal}\n",
+       "(True)\n",
+       "\n",
+       "\n",
+       "Block griddata\n",
+       "--------------------\n",
+       "icelltype\n",
+       "Layer_1{constant 1}\n",
+       "Layer_2{constant 0}\n",
+       "\n",
+       "k\n",
+       "Layer_1{constant 10.0}\n",
+       "Layer_2{constant 5.0}\n",
+       "\n",
+       "k33\n",
+       "{constant 0.1}\n",
+       "\n"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fp.mf6.ModflowGwfnpf(\n",
+    "    model=model,\n",
+    "    save_specific_discharge=True,\n",
+    "    icelltype=[1, 0],\n",
+    "    k=[10.0, 5.0],\n",
+    "    k33=0.1,\n",
+    "    filename=f\"{sim_name}.npf\",\n",
+    "    pname=\"npf\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "package_name = ic\n",
+       "filename = gwf_bcf2ss.ic\n",
+       "package_type = ic\n",
+       "model_or_simulation_package = model\n",
+       "model_name = gwf_bcf2ss\n",
+       "\n",
+       "Block griddata\n",
+       "--------------------\n",
+       "strt\n",
+       "{constant 0.0}\n",
+       "\n"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fp.mf6.ModflowGwfic(\n",
+    "    model=model,\n",
+    "    strt=0.0,\n",
+    "    filename=f\"{sim_name}.ic\",\n",
+    "    pname=\"ic\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "package_name = riv_0\n",
+       "filename = gwf_bcf2ss.riv\n",
+       "package_type = riv\n",
+       "model_or_simulation_package = model\n",
+       "model_name = gwf_bcf2ss\n",
+       "\n",
+       "Block dimensions\n",
+       "--------------------\n",
+       "maxbound\n",
+       "{internal}\n",
+       "(10)\n",
+       "\n",
+       "\n",
+       "Block period\n",
+       "--------------------\n",
+       "stress_period_data\n",
+       "{0: <flopy.mf6.data.mfdataplist.PandasListStorage object at 0x0000028BDE4567D0>}\n"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "def make_riv_iter():\n",
+    "    for i in range(10):\n",
+    "        yield ((1, i, 14), 0.0, 10000.0, -5.0)\n",
+    "\n",
+    "fp.mf6.ModflowGwfriv(\n",
+    "    model=model,\n",
+    "    maxbound=10,\n",
+    "    stress_period_data=list(make_riv_iter()),\n",
+    "    filename=f\"{sim_name}.riv\",\n",
+    "    pname=\"riv_0\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "rec.array([((1, 0, 14), 0., 10000., -5.), ((1, 1, 14), 0., 10000., -5.),\n",
+       "           ((1, 2, 14), 0., 10000., -5.), ((1, 3, 14), 0., 10000., -5.),\n",
+       "           ((1, 4, 14), 0., 10000., -5.), ((1, 5, 14), 0., 10000., -5.),\n",
+       "           ((1, 6, 14), 0., 10000., -5.), ((1, 7, 14), 0., 10000., -5.),\n",
+       "           ((1, 8, 14), 0., 10000., -5.), ((1, 9, 14), 0., 10000., -5.)],\n",
+       "          dtype=[('cellid', 'O'), ('stage', '<f8'), ('cond', '<f8'), ('rbot', '<f8')])"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "model.package_dict['riv_0'].stress_period_data.get_data(key=0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "package_name = wel_0\n",
+       "filename = gwf_bcf2ss.wel\n",
+       "package_type = wel\n",
+       "model_or_simulation_package = model\n",
+       "model_name = gwf_bcf2ss\n",
+       "\n",
+       "Block dimensions\n",
+       "--------------------\n",
+       "maxbound\n",
+       "{internal}\n",
+       "(2)\n",
+       "\n",
+       "\n",
+       "Block period\n",
+       "--------------------\n",
+       "stress_period_data\n",
+       "{1: <flopy.mf6.data.mfdataplist.PandasListStorage object at 0x0000028BDE177B10>}\n"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fp.mf6.ModflowGwfwel(\n",
+    "    model=model,\n",
+    "    maxbound=2,\n",
+    "    stress_period_data={\n",
+    "        1: [\n",
+    "            ((1, 2, 3), -35000.0),\n",
+    "            ((1, 7, 3), -35000.0)\n",
+    "        ]\n",
+    "    },\n",
+    "    filename=f\"{sim_name}.wel\",\n",
+    "    pname=\"wel_0\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "rec.array([((1, 2, 3), -35000.), ((1, 7, 3), -35000.)],\n",
+       "          dtype=[('cellid', 'O'), ('q', '<f8')])"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "model.package_dict['wel_0'].stress_period_data.get_data(key=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "package_name = rcha_0\n",
+       "filename = gwf_bcf2ss.rcha\n",
+       "package_type = rcha\n",
+       "model_or_simulation_package = model\n",
+       "model_name = gwf_bcf2ss\n",
+       "\n",
+       "Block options\n",
+       "--------------------\n",
+       "readasarrays\n",
+       "{internal}\n",
+       "(True)\n",
+       "\n",
+       "\n",
+       "Block period\n",
+       "--------------------\n",
+       "irch\n",
+       "None\n",
+       "recharge\n",
+       "{constant 0.004}\n",
+       "\n",
+       "aux\n",
+       "None\n"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fp.mf6.ModflowGwfrcha(\n",
+    "    model=model,\n",
+    "    recharge=0.004,\n",
+    "    filename=f\"{sim_name}.rcha\",\n",
+    "    pname=\"rcha_0\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "package_name = oc\n",
+       "filename = gwf_bcf2ss.oc\n",
+       "package_type = oc\n",
+       "model_or_simulation_package = model\n",
+       "model_name = gwf_bcf2ss\n",
+       "\n",
+       "Block options\n",
+       "--------------------\n",
+       "budget_filerecord\n",
+       "{internal}\n",
+       "(rec.array([('gwf_bcf2ss.cbc',)],\n",
+       "          dtype=[('budgetfile', 'O')]))\n",
+       "\n",
+       "head_filerecord\n",
+       "{internal}\n",
+       "(rec.array([('gwf_bcf2ss.hds',)],\n",
+       "          dtype=[('headfile', 'O')]))\n",
+       "\n",
+       "\n",
+       "Block period\n",
+       "--------------------\n",
+       "saverecord\n",
+       "{internal}\n",
+       "(rec.array([('head', 'all', None), ('budget', 'all', None)],\n",
+       "          dtype=[('rtype', 'O'), ('ocsetting', 'O'), ('ocsetting_data', 'O')]))\n",
+       "\n",
+       "printrecord\n",
+       "None\n"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fp.mf6.ModflowGwfoc(\n",
+    "    model=model,\n",
+    "    budget_filerecord=[f\"{sim_name}.cbc\"],\n",
+    "    head_filerecord=[f\"{sim_name}.hds\"],\n",
+    "    saverecord=[\n",
+    "        (\"head\", \"all\"),\n",
+    "        (\"budget\", \"all\")\n",
+    "    ],\n",
+    "    filename=f\"{sim_name}.oc\",\n",
+    "    pname=\"oc\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Checking model \"gwf_bcf2ss\"...\n",
+      "\n",
+      "gwf_bcf2ss MODEL DATA VALIDATION SUMMARY:\n",
+      "  No errors or warnings encountered.\n",
+      "\n",
+      "  Checks that passed:\n",
+      "    npf package: zero or negative horizontal hydraulic conductivity values\n",
+      "    npf package: vertical hydraulic conductivity values below checker threshold of 1e-11\n",
+      "    npf package: vertical hydraulic conductivity values above checker threshold of 100000.0\n",
+      "    npf package: horizontal hydraulic conductivity values below checker threshold of 1e-11\n",
+      "    npf package: horizontal hydraulic conductivity values above checker threshold of 100000.0\n",
+      "    riv_0 package: BC indices valid\n",
+      "    riv_0 package: not a number (Nan) entries\n",
+      "    riv_0 package: BC in inactive cells\n",
+      "    wel_0 package: BC indices valid\n",
+      "    wel_0 package: not a number (Nan) entries\n",
+      "    wel_0 package: BC in inactive cells\n",
+      "\n",
+      "Checking for missing simulation packages...\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "[<flopy.utils.check.mf6check at 0x28bde140e10>]"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "sim.check()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "writing simulation...\n",
+      "  writing simulation name file...\n",
+      "  writing simulation tdis package...\n",
+      "  writing solution package ims...\n",
+      "  writing model gwf_bcf2ss...\n",
+      "    writing model name file...\n",
+      "    writing package dis...\n",
+      "    writing package npf...\n",
+      "    writing package ic...\n",
+      "    writing package riv_0...\n",
+      "    writing package wel_0...\n",
+      "    writing package rcha_0...\n",
+      "    writing package oc...\n"
+     ]
+    }
+   ],
+   "source": [
+    "sim.write_simulation()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "FloPy is using the following executable to run the model: ..\\..\\..\\..\\..\\..\\Software\\mf6.4.1\\bin\\mf6.exe\n",
+      "                                   MODFLOW 6\n",
+      "                U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL\n",
+      "                        VERSION 6.4.1 Release 12/09/2022\n",
+      "\n",
+      "   MODFLOW 6 compiled Dec 09 2022 18:18:37 with Intel(R) Fortran Intel(R) 64\n",
+      "   Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0\n",
+      "                             Build 20220726_000000\n",
+      "\n",
+      "This software has been approved for release by the U.S. Geological \n",
+      "Survey (USGS). Although the software has been subjected to rigorous \n",
+      "review, the USGS reserves the right to update the software as needed \n",
+      "pursuant to further analysis and review. No warranty, expressed or \n",
+      "implied, is made by the USGS or the U.S. Government as to the \n",
+      "functionality of the software and related material nor shall the \n",
+      "fact of release constitute any such warranty. Furthermore, the \n",
+      "software is released on condition that neither the USGS nor the U.S. \n",
+      "Government shall be held liable for any damages resulting from its \n",
+      "authorized or unauthorized use. Also refer to the USGS Water \n",
+      "Resources Software User Rights Notice for complete use, copyright, \n",
+      "and distribution information.\n",
+      "\n",
+      " \n",
+      " Run start date and time (yyyy/mm/dd hh:mm:ss): 2024/01/13 16:12:58\n",
+      " \n",
+      " Writing simulation list file: mfsim.lst\n",
+      " Using Simulation name file: mfsim.nam\n",
+      " \n",
+      "    Solving:  Stress period:     1    Time step:     1\n",
+      "    Solving:  Stress period:     2    Time step:     1\n",
+      " \n",
+      " Run end date and time (yyyy/mm/dd hh:mm:ss): 2024/01/13 16:12:58\n",
+      " Elapsed run time:  0.042 Seconds\n",
+      " \n",
+      " Normal termination of simulation.\n"
+     ]
+    }
+   ],
+   "source": [
+    "success, buff = sim.run_simulation()\n",
+    "assert success, \"MODFLOW did not terminate normally!\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim.delete_output_files()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.11.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.dis b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.dis
index 23c2e60..e960612 100644
--- a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.dis
+++ b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.dis
@@ -1,4 +1,4 @@
-# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:17:27.
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:33:36.
 BEGIN options
   LENGTH_UNITS  feet
 END options
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.ic b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.ic
index 8121b0c..c2fea1e 100644
--- a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.ic
+++ b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.ic
@@ -1,4 +1,4 @@
-# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:17:27.
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:33:36.
 BEGIN griddata
   strt
     CONSTANT       0.00000000
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.ims b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.ims
index a920ca2..4f0f12c 100644
--- a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.ims
+++ b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.ims
@@ -1,4 +1,4 @@
-# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:17:27.
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:33:36.
 BEGIN options
 END options
 
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.nam b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.nam
index 108cbc9..86cde66 100644
--- a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.nam
+++ b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.nam
@@ -1,15 +1,15 @@
-# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:17:27.
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:33:36.
 BEGIN options
   SAVE_FLOWS
 END options
 
 BEGIN packages
-  DIS6  gwf_bcf2ss.dis  dis
-  NPF6  gwf_bcf2ss.npf  npf
-  IC6  gwf_bcf2ss.ic  ic
-  RIV6  gwf_bcf2ss.riv  riv_0
-  WEL6  gwf_bcf2ss.wel  wel_0
   RCH6  gwf_bcf2ss.rcha  rcha_0
-  OC6  gwf_bcf2ss  oc
+  WEL6  gwf_bcf2ss.wel  wel_0
+  RIV6  gwf_bcf2ss.riv  riv_0
+  IC6  gwf_bcf2ss.ic  ic
+  NPF6  gwf_bcf2ss.npf  npf
+  DIS6  gwf_bcf2ss.dis  dis
+  OC6  gwf_bcf2ss.oc  oc
 END packages
 
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.npf b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.npf
index 874c434..0a6bd70 100644
--- a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.npf
+++ b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.npf
@@ -1,4 +1,4 @@
-# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:17:27.
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:33:36.
 BEGIN options
   REWET  wetfct       1.00000000  iwetit  1  ihdwet  0
   SAVE_SPECIFIC_DISCHARGE
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.oc
similarity index 70%
rename from notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss
rename to notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.oc
index b7bca86..969bda0 100644
--- a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss
+++ b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.oc
@@ -1,4 +1,4 @@
-# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:17:27.
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:33:36.
 BEGIN options
   BUDGET  FILEOUT  gwf_bcf2ss.cbc
   HEAD  FILEOUT  gwf_bcf2ss.hds
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.rcha b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.rcha
index d47963c..6b7131a 100644
--- a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.rcha
+++ b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.rcha
@@ -1,4 +1,4 @@
-# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:17:27.
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:33:36.
 BEGIN options
   READASARRAYS
 END options
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.riv b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.riv
index 7e139ac..c08016f 100644
--- a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.riv
+++ b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.riv
@@ -1,4 +1,4 @@
-# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:17:27.
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:33:36.
 BEGIN options
 END options
 
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.tdis b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.tdis
index 773d217..9fff2bd 100644
--- a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.tdis
+++ b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.tdis
@@ -1,4 +1,4 @@
-# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:17:27.
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:33:36.
 BEGIN options
   TIME_UNITS  days
 END options
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.wel b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.wel
index 7fbf99f..d1405c7 100644
--- a/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.wel
+++ b/notebooks/implementing/gwf_bcf2ss/model_p1a/gwf_bcf2ss.wel
@@ -1,4 +1,4 @@
-# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:17:27.
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:33:36.
 BEGIN options
 END options
 
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p1a/mfsim.nam b/notebooks/implementing/gwf_bcf2ss/model_p1a/mfsim.nam
index 1288453..d78b418 100644
--- a/notebooks/implementing/gwf_bcf2ss/model_p1a/mfsim.nam
+++ b/notebooks/implementing/gwf_bcf2ss/model_p1a/mfsim.nam
@@ -1,4 +1,4 @@
-# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:17:27.
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 15:33:36.
 BEGIN options
 END options
 
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.dis b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.dis
new file mode 100644
index 0000000..fb3fdfa
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.dis
@@ -0,0 +1,23 @@
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+BEGIN options
+  LENGTH_UNITS  feet
+END options
+
+BEGIN dimensions
+  NLAY  2
+  NROW  10
+  NCOL  15
+END dimensions
+
+BEGIN griddata
+  delr
+    CONSTANT     500.00000000
+  delc
+    CONSTANT     500.00000000
+  top
+    CONSTANT     150.00000000
+  botm  LAYERED
+    CONSTANT      50.00000000
+    CONSTANT     -50.00000000
+END griddata
+
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.dis.grb b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.dis.grb
new file mode 100644
index 0000000000000000000000000000000000000000..b659b72ec73c21cd9fdd5e361100f26cfbd99720
GIT binary patch
literal 16048
zcmeI&36$3J8o=?977djoO48fDt3g7#nwpWxG{ZEhh)P02LPDry3keA+OSbIUvM(iD
zLXsq*k|YTU$=3aT>;Ig4#JT34GwyNEeLLrKpJ#saAMf;?ci#Rp=XNg2%Sq+rbxW1>
zPv&8{McwiW3sPAn{)mF&p2exGy-Rw1DwUa^yZK{Uw=UraWfpYG&Cl<a$}1?&?VMYb
zD#*#}nrfY@pW3s;|A{|+q|Ac+>|VcrggyTB1OI~yG7E|dd;I<}vP#;aBo9Ta!u&sY
zj5a0t;1XVwS<of>4<4dz>3ubLmh>zv%Iln0kjg3S-Z4M-Ute=S$JV>I{>`5fe%@PN
zmf7pS@EF+zor?;K|H{8cDl;cHzvwrQkd?|R%`Z6Ne2h+8K1OLiaea*9!b7%vh&HAD
z7#$0XyKeCyrT8+k=_!_(x5c+d7Ozw#ens8#S^CFsXKwi|$d|lSzXAPv<mD80&E|XX
zH!r;YWlHahl9}ggdvUKra{szV$!xJzMn;*8jEqcF!R_VwpMUqK!-x)>-u=HmXB@rx
z?~$JA4nCpN&x7>y;P3oA$jICL^-0fk2cM_u=Rx{;@b`Wm%-HnlW#pyj=C4z_L%IV6
z(%%o8o(6xT4*Uq5Q4aAl@Q&;?u^Z~630h(gw8j3&#)0U9LUhLw=!fHQ5>CTe7=a6M
z8Lq;0xDmHvD(=T4coNU!CCtG*EX0Rcj#c;y-{D7WL^(eBZBZ3Fqc-ZJDO#Zo+94aA
z(G^ALfg{l$12GswF&v|CF~(v%CgNt?fqO6w)A2NB;uXAs`B;P{Sb^18i}m;gW%xxy
z1yn*c)I=RLL~~@JE!v|K@=%EGI2?U20D~|D!*DJ}V+_V&0w!S!rs94)jK}dTW??qw
zVgVLo8CGHq)?ouSqTJRzCr}yHPy@A54~@|RS=bxx&;hyVf<kmd5A;Sq9E(9X1w(N*
zM&JUB!4(*f>oFNqFctUX5j=_K@e<}>9v0$5EXOK*g>SI|8&R$-&kt0^&Zv$0Xo^;7
zgLcS9XXN8h^uUqmkAWDBp%{))xENzG9usjhZpYo2hUs`3Gw}-Ez<eyi60E>#ti^i#
zjIH=4tbj_WhMK5@hG>o~v_*S#LLLfHj9%!20T_fK7>08(8e=dH6EF!=FctUVAv}&}
zF$=RX7Ync$%diq_unrrr5oP&ZSw&Pqb<{##G)7D8iG6VZa&a&YMGy2wUktzr7>v_#
z7S6+HT!OK<8rNYmrr-|n?X>ClcH1+3J>luP>5MOVMpax7ch5<iu>&S1yXSoR8P#w@
zvU?5x*$iL%CMCPqR6uo1PIj-^20OvmFn6!1h@Ej$vU~39X${<*>|V1id~LlY*}Y~v
z?1CxD?lqNA3%4e_*Hp%?xGmYerV47~_GI@OU$^goZR})^I=D00ZQ}cEl7g?n?$ORJ
zb#WKmqn~Z+;cmD`L;KXnJ#deXHfn%-;T|pR)DZW<J$l-z5&ZYSJ(}99F&=<>bhTL%
zOoMy0wOdm>2>0k~yJmO@?$OwO&G9hYqq7}bz<-<EqjgL6M`2Strm#DvCwug^VJkca
z_vmfI)_5H5(c6YucmnRx+lG7KNw`OE8}5mx;2yngxEG#=d-S$p8$1K|=xxKj@hse<
zw+-9kIk-n}8}5VW;U2whxG!eFJ$l=4Kg@)C^tNL=ya4y;y+8Xb*wn^0Y>yX{J;n5W
zHv9)(O7`e&!vpX#+@rS*JKz<#M{gTu<5jpvZyR>RY`8~n8+O8LaF5<L%)uPEM{gVE
z;&r%3ZyO$nH{c$<ZP*!a!aaK1@F2{Ed-S$p9^Qg`^tNLc%!7OMKA8P&*wl_T%*Xs>
zkKQ)yig(~1y*scMU;*5tw+##NF5IKH4G+OWxJPdr9*Xzi9=&Z?ghg<V-Zt!p_u(GB
zZCH%OaF5<L?2ZrM9=&aN7(Rr1^tNFSd<6ICZNr{e0{7@`$6i<p_vn2%`!d+n#x^_x
zA18bCwqb89hkNw4;gR?R?j7L!ZFm$`BzyF>VIO=7_vmfIzE}zO=xxJ(_zdpR+lEJD
z72Kn@4g2GBxJPdr9)s0zkKQ&MfG^-4y={0b*1$b_+weGi3HRu2$K&x8+@tqE_OD@6
z8{6;%tWEalZNowM2JX?@h9_bj+@rS*Pr|owkKTRRPsVp}kKQ&MjPKzdy={03*26t|
z+i(befP3_|;i=dF_vmfI)9@qQqqhxD$4_vN-ZmVHpWz<8ZFmNLfqV3}<1qXR_vn2l
z`$pK*#x^_)-b}bhZyTQNG3V~l+lIqYhI9AmZNqc0HRtZp+lJ?&Ea&ch;rnfP9{d^4
zy&rtP4M(6n=kC$lh9gmdbNA?N!%^6VbNA?N!}C#*bNA?N!_mm(+&y~R@B;YrrhD|Z
z;f2_abNA?N$BR&jbNA?dF?(gsY-(d0UV<u|xkqmsj=}bvyGL&uUJ8Hab&uXQybL>V
z?jF5ucsZ(Z?jF5uI2Jo{?jF5IvR{GfoV!QwQS9Te6X)*H+lE)dpTFIsw+*jC4bI)8
zw+*jGP0rn;w++W*7tY<Iw+*jBEzaGew;d;7SI*s|_qFW)9>AtHw&8BP?Tps;tHYTy
zI@>PAnKK&OtuANI=xeijoH?Vdz3Ox3jIOq7z?pL-e4m{fa^@Tb-)EynoH?VVeHwG-
zjE=Tx!kIG~+Qr{>IHR9UnsMfgb~b6wnKQci|9!ZQw;7v$Ub-`xWtg5%#MESN&Ac8r
z;I3qrWqQpd+?~vFOs|=Ydy-k6=`}av-egu_dd*F^FPYmgz2;`zpUjF(uek*eBr}ug
zHB&GxncFgN!L4{OncFeF<~BT(%t}nJxg8HDvoh0b?!Y6-titq~M`4%kVH0~yPi9r7
zO?<yi9!us9%xGtq$KkAo=x3WJ;M@_>&^}MXSsl^QMo+=H6QZS^o`!R0L{D2i17{6H
zQ+quNXHDD!ud&&4$=rn*ZSD3voV5^rZ8rnXT@j7#Hxtg<h|YF=0nXhJt!J^@vkvTN
z!xxj8Vn%Npz657oL~k3u3}-z=ZyUY>XMIF(8@>u>14M5d&W5ugqPGoSgR>E$w+-jO
z*%;B=hOfif1ku}uZ@}3U(R&8_n{YNm^tRz#IGZDS+wd(oTOfMdaUPs45xw7Lx99G#
zvkm7bvlTOX+wdJYTO)egZ~>fIh~74Q7tTEpy=}M<&OH&mZTKFXdm(z;a1oqs5WQ{q
zKAd|adfRX@oNW=kZTJD4`yhJT@IyHF#SD0j4L?fee$41?!zFOGL-e-eQaJZV^j^kp
z&-Sph4L?riKbX<mhRfkR0MXlqpTOAx(c6YA;LJw!w&AC6c0}~H;Yv6=A$r^JGdOb)
zy=}M(&Rj%q8-5Pwfr#EVTn%SuL~k2@0p~%8-ZoqVXC9*WQuZ(5?1JcR$FJZ#7}5J{
zc6;W-&Nf_|%&yGnZNqQiEI{<O;W{`A5xs5rEu4oSdfV_jI1feiw&C}179o1ua6O#e
z5WQ{q1DwT(-ZtC-XLm$z8~zCAVTj%~{0Yt;h~75*8P1-F-ZuON&R$pwud(B=$vm7H
zy*IMk^9b14hTaf3dn0<=(EEJnk%-<l^giEt6r#5cz0Y^{LG-p^S<an(5xs3#j&o-}
zL~k3G=iGTTqPGpb&v*7m^tPe*`Oaeyy=~}yzH<Pgw++3|cOHxAZA0(#oyQ@1+tB-b
z=kbW%cJw~qIS|p?`+R$zz_7Coz0Y?J;_c{dL+|sQCvqOWZRmZz^CZrrw++3|cb?37
z^tPe*`Od+dM{gT?pYJ?{^XP3u@AI8QIFH^o^giEtD(BJLhTi8pPvbm#+tB-b=johB
zZyS1_?;Of`^tPe*`OY&qkKQ))KHoWv^XP3y@AI8!avr_C&$s7UoY~oi-s?Nh_K4~G
zZRq{Jb2#VG+lDDP&*40J+psR2=W-ssZCDS^^Ei**Hmncl2+pIo4I98YlJn?o!-jB<
z;yilWuo0Z+a~{2I*ci^yoJVgPHi7d3&ZD;to5Fb^=h54S&EUL<^XP5I=5SujdGv0<
z+xEPKGdtQ)V|Y7S+fSEr7M*RU%Q%b1cGKmYMPHj~EN9WyUb=#_=xQsC<1CumNmp_f
zJ#D0`IE$9{(bb$qN84yTXVK6ux`wmpXA@1}EZW&b*K!u!crzm{-5~uONIwVC&w=!F
S;Q#eGz(0fOKj{LcUf`bzQ*E{Y

literal 0
HcmV?d00001

diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.ic b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.ic
new file mode 100644
index 0000000..87c4323
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.ic
@@ -0,0 +1,6 @@
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+BEGIN griddata
+  strt
+    CONSTANT       0.00000000
+END griddata
+
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.ims b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.ims
new file mode 100644
index 0000000..3739ebe
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.ims
@@ -0,0 +1,17 @@
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+BEGIN options
+END options
+
+BEGIN nonlinear
+  OUTER_DVCLOSE  1.00000000E-06
+  OUTER_MAXIMUM  500
+END nonlinear
+
+BEGIN linear
+  INNER_MAXIMUM  100
+  INNER_DVCLOSE  1.00000000E-06
+  inner_rclose       0.00100000  strict
+  LINEAR_ACCELERATION  bicgstab
+  RELAXATION_FACTOR       0.97000000
+END linear
+
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.lst b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.lst
new file mode 100644
index 0000000..161afe0
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.lst
@@ -0,0 +1,377 @@
+                                   MODFLOW 6
+                U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
+                          GROUNDWATER FLOW MODEL (GWF)
+                        VERSION 6.4.1 Release 12/09/2022
+
+   MODFLOW 6 compiled Dec 09 2022 18:18:37 with Intel(R) Fortran Intel(R) 64
+   Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0
+                             Build 20220726_000000
+
+This software has been approved for release by the U.S. Geological 
+Survey (USGS). Although the software has been subjected to rigorous 
+review, the USGS reserves the right to update the software as needed 
+pursuant to further analysis and review. No warranty, expressed or 
+implied, is made by the USGS or the U.S. Government as to the 
+functionality of the software and related material nor shall the 
+fact of release constitute any such warranty. Furthermore, the 
+software is released on condition that neither the USGS nor the U.S. 
+Government shall be held liable for any damages resulting from its 
+authorized or unauthorized use. Also refer to the USGS Water 
+Resources Software User Rights Notice for complete use, copyright, 
+and distribution information.
+
+
+As a work of the United States Government, this USGS product is 
+in the public domain within the United States. You can copy, 
+modify, distribute, and perform the work, even for commercial 
+purposes, all without asking permission. Additionally, USGS 
+waives copyright and related rights in the work worldwide 
+through CC0 1.0 Universal Public Domain Dedication 
+(https://creativecommons.org/publicdomain/zero/1.0/).
+
+The following GNU Lesser General Public License (LGPL) libraries
+are used in this USGS product:
+
+    SPARSKIT version 2.0
+      ilut, luson, and qsplit 
+      (https://www-users.cse.umn.edu/~saad/software/SPARSKIT/)
+
+    RCM - Reverse Cuthill McKee Ordering
+      (https://people.math.sc.edu/Burkardt/f_src/rcm/rcm.html)
+
+    BLAS - Basic Linear Algebra Subprograms Level 1
+      (https://people.math.sc.edu/Burkardt/f_src/blas1_d/blas1_d.html)
+
+    SPARSEKIT - Sparse Matrix Utility Package
+      amux, dperm, dvperm, rperm, and cperm
+      (https://people.sc.fsu.edu/~jburkardt/f77_src/sparsekit/sparsekit.html)
+
+The following BSD-3 License libraries are used in this USGS product:
+
+    Modern Fortran DAG Library
+      Copyright (c) 2018, Jacob Williams
+      All rights reserved.
+      (https://github.com/jacobwilliams/daglib)
+
+MODFLOW 6 compiler options: /Isrc\libmf6core.a.p /Isrc /I..\src
+/Isrc\libmf6_external.a.p /MTd /nologo /warn:general /warn:truncated_source
+/stand:f08 /Zi /traceback /O2 /fpe:0 /heap-arrays:0 /traceback /fpp
+/Qdiag-disable:7416 /Qdiag-disable:7025 /Qdiag-disable:5268
+/Fdsrc\libmf6core.a.p\Utilities_compilerversion.F90.pdb
+/module:src\libmf6core.a.p
+/Fosrc/libmf6core.a.p/Utilities_compilerversion.F90.obj /c
+
+System command used to initiate simulation:
+C:\Users\buckl\Documents\Software\mf6.4.1\bin\mf6.exe
+
+MODFLOW was compiled using uniform precision.
+
+Real Variables
+  KIND: 8
+  TINY (smallest non-zero value):    2.225074-308
+  HUGE (largest value):    1.797693+308
+  PRECISION: 15
+  SIZE IN BITS: 64
+
+Integer Variables
+  KIND: 4
+  HUGE (largest value): 2147483647
+  SIZE IN BITS: 32
+
+Long Integer Variables
+  KIND: 8
+  HUGE (largest value): 9223372036854775807
+  SIZE IN BITS: 64
+
+Logical Variables
+  KIND: 4
+  SIZE IN BITS: 32
+ 
+
+ OPENED gwf_bcf2ss.dis
+ FILE TYPE:DIS6                   UNIT 1005   STATUS:OLD                 
+ FORMAT:FORMATTED              ACCESS:SEQUENTIAL          
+ ACTION:READ                
+
+
+ OPENED gwf_bcf2ss.npf
+ FILE TYPE:NPF6                   UNIT 1006   STATUS:OLD                 
+ FORMAT:FORMATTED              ACCESS:SEQUENTIAL          
+ ACTION:READ                
+
+
+ OPENED gwf_bcf2ss.ic
+ FILE TYPE:IC6                    UNIT 1007   STATUS:OLD                 
+ FORMAT:FORMATTED              ACCESS:SEQUENTIAL          
+ ACTION:READ                
+
+
+ OPENED gwf_bcf2ss.riv
+ FILE TYPE:RIV6                   UNIT 1008   STATUS:OLD                 
+ FORMAT:FORMATTED              ACCESS:SEQUENTIAL          
+ ACTION:READ                
+
+
+ OPENED gwf_bcf2ss.wel
+ FILE TYPE:WEL6                   UNIT 1009   STATUS:OLD                 
+ FORMAT:FORMATTED              ACCESS:SEQUENTIAL          
+ ACTION:READ                
+
+
+ OPENED gwf_bcf2ss.rcha
+ FILE TYPE:RCH6                   UNIT 1010   STATUS:OLD                 
+ FORMAT:FORMATTED              ACCESS:SEQUENTIAL          
+ ACTION:READ                
+
+
+ OPENED gwf_bcf2ss.oc
+ FILE TYPE:OC6                    UNIT 1011   STATUS:OLD                 
+ FORMAT:FORMATTED              ACCESS:SEQUENTIAL          
+ ACTION:READ                
+
+ NAMEFILE OPTIONS:
+    FLOWS WILL BE SAVED TO BUDGET FILE SPECIFIED IN OUTPUT CONTROL
+    NEWTON-RAPHSON method enabled for the model.
+
+ DIS -- STRUCTURED GRID DISCRETIZATION PACKAGE, VERSION 2 : 3/27/2014 - INPUT READ FROM UNIT 1005
+
+ Loading input for GWF_BCF2SS/DIS
+ # File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+   NLAY = 2
+   NROW = 10
+   NCOL = 15
+   DELR = 500.0000000000000
+   DELC = 500.0000000000000
+   TOP = 150.0000000000000
+   BOTM = variable 3D double precision array ranging from -50.00000000000000 to 50.00000000000000
+ Loading input complete...
+
+
+ NPF -- NODE PROPERTY FLOW PACKAGE, VERSION 1, 3/30/2015 INPUT READ FROM UNIT 1006
+
+ Loading input for GWF_BCF2SS/NPF
+ # File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+   ISAVSPDIS = 1
+   ICELLTYPE = variable 1D integer array ranging from 0 to 1
+   K = variable 1D double precision array ranging from 5.000000000000000 to 10.00000000000000
+   K33 = .1000000000000000
+ Loading input complete...
+
+ Setting Discretization Options
+    Model length unit [0=UND, 1=FEET, 2=METERS, 3=CENTIMETERS] set as 1
+ End Setting Discretization Options
+
+ Setting Discretization Dimensions
+    NLAY = 2
+    NROW = 10
+    NCOL = 15
+ End Setting Discretization Dimensions
+
+ Setting Discretization Griddata
+    DELR set from input file
+    DELC set from input file
+    TOP set from input file
+    BOTM set from input file
+ End Setting Discretization Griddata
+
+ Setting NPF Options
+    Specific discharge will be calculated at cell centers and written to DATA-SPDIS in budget file when requested.
+ End Setting NPF Options
+
+ Setting NPF Griddata
+    ICELLTYPE set from input file
+    K set from input file
+    K33 set from input file
+    K22 not provided.  Setting K22 = K.
+ End Setting NPF Griddata
+
+
+ WEL   -- WEL PACKAGE, VERSION 8, 2/22/2014 INPUT READ FROM UNIT 1009
+ # File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:16.
+
+ PROCESSING WEL OPTIONS
+ END OF WEL OPTIONS
+
+ PROCESSING WEL DIMENSIONS
+    MAXBOUND =       2
+ END OF WEL DIMENSIONS
+
+ RIV   -- RIV PACKAGE, VERSION 8, 2/22/2014 INPUT READ FROM UNIT 1008
+ # File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+
+ PROCESSING RIV OPTIONS
+ END OF RIV OPTIONS
+
+ PROCESSING RIV DIMENSIONS
+    MAXBOUND =      10
+ END OF RIV DIMENSIONS
+
+ RCH   -- RCH PACKAGE, VERSION 8, 2/22/2014 INPUT READ FROM UNIT 1010
+ # File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:16.
+
+ PROCESSING RCH OPTIONS
+    RECHARGE INPUT WILL BE READ AS ARRAY(S).
+ END OF RCHA OPTIONS
+
+ IC -- INITIAL CONDITIONS PACKAGE, VERSION 8, 3/28/2015 INPUT READ FROM UNIT 1007
+ # File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+ PROCESSING GRIDDATA
+
+     INITIAL HEAD =   0.000000    
+ END PROCESSING GRIDDATA
+    BINARY GRID INFORMATION WILL BE WRITTEN TO:
+      UNIT NUMBER: 1013
+      FILE NAME: c:\Users\buckl\Documents\Github\flopy_mf6_work\notebooks\implementing\gwf_bcf2ss\model_p2a\gwf_bcf2ss.dis.grb
+
+ OPENED c:\Users\buckl\Documents\Github\flopy_mf6_work\notebooks\implementing\gwf_bcf2ss\model_p2a\gwf_bcf2ss.dis.grb
+ FILE TYPE:DATA(BINARY)   UNIT 1013   STATUS:REPLACE             
+ FORMAT:UNFORMATTED            ACCESS:STREAM              
+ ACTION:READWRITE           
+
+      THE LAST TIME STEP WILL BE PRINTED
+      THE LAST TIME STEP WILL BE PRINTED
+ # File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:16.
+
+ PROCESSING OC OPTIONS
+
+    BUDGET INFORMATION WILL BE WRITTEN TO:
+      UNIT NUMBER: 1014
+      FILE NAME: gwf_bcf2ss.cbc
+
+ OPENED gwf_bcf2ss.cbc
+ FILE TYPE:DATA(BINARY)   UNIT 1014   STATUS:REPLACE             
+ FORMAT:UNFORMATTED            ACCESS:STREAM              
+ ACTION:READWRITE           
+
+    HEAD INFORMATION WILL BE WRITTEN TO:
+      UNIT NUMBER: 1015
+      FILE NAME: gwf_bcf2ss.hds
+
+ OPENED gwf_bcf2ss.hds
+ FILE TYPE:DATA(BINARY)   UNIT 1015   STATUS:REPLACE             
+ FORMAT:UNFORMATTED            ACCESS:STREAM              
+ ACTION:READWRITE           
+
+ END OF OC OPTIONS
+
+start timestep kper="1" kstp="1" mode="normal"
+
+
+ BEGIN READING OUTPUT CONTROL FOR STRESS PERIOD 1
+      ALL TIME STEPS WILL BE SAVED
+      ALL TIME STEPS WILL BE SAVED
+
+ END READING OUTPUT CONTROL FOR STRESS PERIOD 1
+
+ REUSING WELS FROM LAST STRESS PERIOD
+
+                 RECHARGE =  0.4000000E-02
+ UBDSV1 SAVING     FLOW-JA-FACE ON UNIT   1014 AT TIME STEP      1, STRESS PERIOD      1
+ UBDSV06 SAVING       DATA-SPDIS IN MODEL GWF_BCF2SS       PACKAGE NPF             CONNECTED TO MODEL GWF_BCF2SS       PACKAGE NPF              ON UNIT   1014 AT TIME STEP      1, STRESS PERIOD      1
+ UBDSV06 SAVING              WEL IN MODEL GWF_BCF2SS       PACKAGE GWF_BCF2SS      CONNECTED TO MODEL GWF_BCF2SS       PACKAGE WEL_0            ON UNIT   1014 AT TIME STEP      1, STRESS PERIOD      1
+ UBDSV06 SAVING              RIV IN MODEL GWF_BCF2SS       PACKAGE GWF_BCF2SS      CONNECTED TO MODEL GWF_BCF2SS       PACKAGE RIV_0            ON UNIT   1014 AT TIME STEP      1, STRESS PERIOD      1
+ UBDSV06 SAVING             RCHA IN MODEL GWF_BCF2SS       PACKAGE GWF_BCF2SS      CONNECTED TO MODEL GWF_BCF2SS       PACKAGE RCHA_0           ON UNIT   1014 AT TIME STEP      1, STRESS PERIOD      1
+
+ HEAD WILL BE SAVED ON UNIT 1015 AT END OF TIME STEP    1, STRESS PERIOD    1
+
+
+  VOLUME BUDGET FOR ENTIRE MODEL AT END OF TIME STEP    1, STRESS PERIOD   1
+  ---------------------------------------------------------------------------------------------------
+
+     CUMULATIVE VOLUME      L**3       RATES FOR THIS TIME STEP      L**3/T          PACKAGE NAME    
+     ------------------                 ------------------------                     ----------------
+
+           IN:                                      IN:
+           ---                                      ---
+                 WEL =           0.0000                   WEL =           0.0000     WEL_0                            
+                 RIV =           0.0000                   RIV =           0.0000     RIV_0                            
+                RCHA =      150000.0000                  RCHA =      150000.0000     RCHA_0                           
+
+            TOTAL IN =      150000.0000              TOTAL IN =      150000.0000
+
+          OUT:                                     OUT:
+          ----                                     ----
+                 WEL =           0.0000                   WEL =           0.0000     WEL_0                            
+                 RIV =      149999.9978                   RIV =      149999.9978     RIV_0                            
+                RCHA =           0.0000                  RCHA =           0.0000     RCHA_0                           
+
+           TOTAL OUT =      149999.9978             TOTAL OUT =      149999.9978
+
+            IN - OUT =       2.1948E-03              IN - OUT =       2.1948E-03
+
+ PERCENT DISCREPANCY =           0.00     PERCENT DISCREPANCY =           0.00
+
+
+
+
+         TIME SUMMARY AT END OF TIME STEP    1 IN STRESS PERIOD    1
+                    SECONDS     MINUTES      HOURS       DAYS        YEARS
+                    -----------------------------------------------------------
+   TIME STEP LENGTH  86400.      1440.0      24.000      1.0000     2.73785E-03
+ STRESS PERIOD TIME  86400.      1440.0      24.000      1.0000     2.73785E-03
+         TOTAL TIME  86400.      1440.0      24.000      1.0000     2.73785E-03
+
+
+end timestep
+
+
+start timestep kper="2" kstp="1" mode="normal"
+
+
+ END OF FILE DETECTED IN OUTPUT CONTROL.
+ CURRENT OUTPUT CONTROL SETTINGS WILL BE 
+ REPEATED UNTIL THE END OF THE SIMULATION.
+
+ OUTPUT CONTROL FOR STRESS PERIOD 2 IS REPEATED USING SETTINGS FROM A PREVIOUS STRESS PERIOD.
+
+ REUSING RIVS FROM LAST STRESS PERIOD
+
+ REUSING RCHS FROM LAST STRESS PERIOD
+ UBDSV1 SAVING     FLOW-JA-FACE ON UNIT   1014 AT TIME STEP      1, STRESS PERIOD      2
+ UBDSV06 SAVING       DATA-SPDIS IN MODEL GWF_BCF2SS       PACKAGE NPF             CONNECTED TO MODEL GWF_BCF2SS       PACKAGE NPF              ON UNIT   1014 AT TIME STEP      1, STRESS PERIOD      2
+ UBDSV06 SAVING              WEL IN MODEL GWF_BCF2SS       PACKAGE GWF_BCF2SS      CONNECTED TO MODEL GWF_BCF2SS       PACKAGE WEL_0            ON UNIT   1014 AT TIME STEP      1, STRESS PERIOD      2
+ UBDSV06 SAVING              RIV IN MODEL GWF_BCF2SS       PACKAGE GWF_BCF2SS      CONNECTED TO MODEL GWF_BCF2SS       PACKAGE RIV_0            ON UNIT   1014 AT TIME STEP      1, STRESS PERIOD      2
+ UBDSV06 SAVING             RCHA IN MODEL GWF_BCF2SS       PACKAGE GWF_BCF2SS      CONNECTED TO MODEL GWF_BCF2SS       PACKAGE RCHA_0           ON UNIT   1014 AT TIME STEP      1, STRESS PERIOD      2
+
+ HEAD WILL BE SAVED ON UNIT 1015 AT END OF TIME STEP    1, STRESS PERIOD    2
+
+
+  VOLUME BUDGET FOR ENTIRE MODEL AT END OF TIME STEP    1, STRESS PERIOD   2
+  ---------------------------------------------------------------------------------------------------
+
+     CUMULATIVE VOLUME      L**3       RATES FOR THIS TIME STEP      L**3/T          PACKAGE NAME    
+     ------------------                 ------------------------                     ----------------
+
+           IN:                                      IN:
+           ---                                      ---
+                 WEL =           0.0000                   WEL =           0.0000     WEL_0                            
+                 RIV =           0.0000                   RIV =           0.0000     RIV_0                            
+                RCHA =      300000.0000                  RCHA =      150000.0000     RCHA_0                           
+
+            TOTAL IN =      300000.0000              TOTAL IN =      150000.0000
+
+          OUT:                                     OUT:
+          ----                                     ----
+                 WEL =       70000.0000                   WEL =       70000.0000     WEL_0                            
+                 RIV =      229999.9978                   RIV =       80000.0000     RIV_0                            
+                RCHA =           0.0000                  RCHA =           0.0000     RCHA_0                           
+
+           TOTAL OUT =      299999.9978             TOTAL OUT =      150000.0000
+
+            IN - OUT =       2.1556E-03              IN - OUT =      -3.9223E-05
+
+ PERCENT DISCREPANCY =           0.00     PERCENT DISCREPANCY =          -0.00
+
+
+
+
+         TIME SUMMARY AT END OF TIME STEP    1 IN STRESS PERIOD    2
+                    SECONDS     MINUTES      HOURS       DAYS        YEARS
+                    -----------------------------------------------------------
+   TIME STEP LENGTH  86400.      1440.0      24.000      1.0000     2.73785E-03
+ STRESS PERIOD TIME  86400.      1440.0      24.000      1.0000     2.73785E-03
+         TOTAL TIME 1.72800E+05  2880.0      48.000      2.0000     5.47570E-03
+
+
+end timestep
+
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.nam b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.nam
new file mode 100644
index 0000000..5f134e3
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.nam
@@ -0,0 +1,16 @@
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+BEGIN options
+  SAVE_FLOWS
+  NEWTON
+END options
+
+BEGIN packages
+  DIS6  gwf_bcf2ss.dis  dis
+  NPF6  gwf_bcf2ss.npf  npf
+  IC6  gwf_bcf2ss.ic  ic
+  RIV6  gwf_bcf2ss.riv  riv_0
+  WEL6  gwf_bcf2ss.wel  wel_0
+  RCH6  gwf_bcf2ss.rcha  rcha_0
+  OC6  gwf_bcf2ss.oc  oc
+END packages
+
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.npf b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.npf
new file mode 100644
index 0000000..b3eb4bf
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.npf
@@ -0,0 +1,16 @@
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+BEGIN options
+  SAVE_SPECIFIC_DISCHARGE
+END options
+
+BEGIN griddata
+  icelltype  LAYERED
+    CONSTANT  1
+    CONSTANT  0
+  k  LAYERED
+    CONSTANT      10.00000000
+    CONSTANT       5.00000000
+  k33
+    CONSTANT       0.10000000
+END griddata
+
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.oc b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.oc
new file mode 100644
index 0000000..a974895
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.oc
@@ -0,0 +1,11 @@
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:16.
+BEGIN options
+  BUDGET  FILEOUT  gwf_bcf2ss.cbc
+  HEAD  FILEOUT  gwf_bcf2ss.hds
+END options
+
+BEGIN period  1
+  SAVE  head  all
+  SAVE  budget  all
+END period  1
+
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.rcha b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.rcha
new file mode 100644
index 0000000..cbb0911
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.rcha
@@ -0,0 +1,10 @@
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:16.
+BEGIN options
+  READASARRAYS
+END options
+
+BEGIN period  1
+  recharge
+    CONSTANT       0.00400000
+END period  1
+
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.riv b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.riv
new file mode 100644
index 0000000..ccf6283
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.riv
@@ -0,0 +1,21 @@
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+BEGIN options
+END options
+
+BEGIN dimensions
+  MAXBOUND  10
+END dimensions
+
+BEGIN period  1
+  2 1 15 0.00000000E+00 1.00000000E+04 -5.00000000E+00
+  2 2 15 0.00000000E+00 1.00000000E+04 -5.00000000E+00
+  2 3 15 0.00000000E+00 1.00000000E+04 -5.00000000E+00
+  2 4 15 0.00000000E+00 1.00000000E+04 -5.00000000E+00
+  2 5 15 0.00000000E+00 1.00000000E+04 -5.00000000E+00
+  2 6 15 0.00000000E+00 1.00000000E+04 -5.00000000E+00
+  2 7 15 0.00000000E+00 1.00000000E+04 -5.00000000E+00
+  2 8 15 0.00000000E+00 1.00000000E+04 -5.00000000E+00
+  2 9 15 0.00000000E+00 1.00000000E+04 -5.00000000E+00
+  2 10 15 0.00000000E+00 1.00000000E+04 -5.00000000E+00
+END period  1
+
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.tdis b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.tdis
new file mode 100644
index 0000000..8db14ba
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.tdis
@@ -0,0 +1,14 @@
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+BEGIN options
+  TIME_UNITS  days
+END options
+
+BEGIN dimensions
+  NPER  2
+END dimensions
+
+BEGIN perioddata
+       1.00000000  1       1.00000000
+       1.00000000  1       1.00000000
+END perioddata
+
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.wel b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.wel
new file mode 100644
index 0000000..697cefd
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/gwf_bcf2ss.wel
@@ -0,0 +1,13 @@
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:16.
+BEGIN options
+END options
+
+BEGIN dimensions
+  MAXBOUND  2
+END dimensions
+
+BEGIN period  2
+  2 3 4 -3.50000000E+04
+  2 8 4 -3.50000000E+04
+END period  2
+
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/mfsim.lst b/notebooks/implementing/gwf_bcf2ss/model_p2a/mfsim.lst
new file mode 100644
index 0000000..15120f5
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/mfsim.lst
@@ -0,0 +1,234 @@
+                                   MODFLOW 6
+                U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
+                        VERSION 6.4.1 Release 12/09/2022
+
+   MODFLOW 6 compiled Dec 09 2022 18:18:37 with Intel(R) Fortran Intel(R) 64
+   Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0
+                             Build 20220726_000000
+
+This software has been approved for release by the U.S. Geological 
+Survey (USGS). Although the software has been subjected to rigorous 
+review, the USGS reserves the right to update the software as needed 
+pursuant to further analysis and review. No warranty, expressed or 
+implied, is made by the USGS or the U.S. Government as to the 
+functionality of the software and related material nor shall the 
+fact of release constitute any such warranty. Furthermore, the 
+software is released on condition that neither the USGS nor the U.S. 
+Government shall be held liable for any damages resulting from its 
+authorized or unauthorized use. Also refer to the USGS Water 
+Resources Software User Rights Notice for complete use, copyright, 
+and distribution information.
+
+
+As a work of the United States Government, this USGS product is 
+in the public domain within the United States. You can copy, 
+modify, distribute, and perform the work, even for commercial 
+purposes, all without asking permission. Additionally, USGS 
+waives copyright and related rights in the work worldwide 
+through CC0 1.0 Universal Public Domain Dedication 
+(https://creativecommons.org/publicdomain/zero/1.0/).
+
+The following GNU Lesser General Public License (LGPL) libraries
+are used in this USGS product:
+
+    SPARSKIT version 2.0
+      ilut, luson, and qsplit 
+      (https://www-users.cse.umn.edu/~saad/software/SPARSKIT/)
+
+    RCM - Reverse Cuthill McKee Ordering
+      (https://people.math.sc.edu/Burkardt/f_src/rcm/rcm.html)
+
+    BLAS - Basic Linear Algebra Subprograms Level 1
+      (https://people.math.sc.edu/Burkardt/f_src/blas1_d/blas1_d.html)
+
+    SPARSEKIT - Sparse Matrix Utility Package
+      amux, dperm, dvperm, rperm, and cperm
+      (https://people.sc.fsu.edu/~jburkardt/f77_src/sparsekit/sparsekit.html)
+
+The following BSD-3 License libraries are used in this USGS product:
+
+    Modern Fortran DAG Library
+      Copyright (c) 2018, Jacob Williams
+      All rights reserved.
+      (https://github.com/jacobwilliams/daglib)
+
+MODFLOW 6 compiler options: /Isrc\libmf6core.a.p /Isrc /I..\src
+/Isrc\libmf6_external.a.p /MTd /nologo /warn:general /warn:truncated_source
+/stand:f08 /Zi /traceback /O2 /fpe:0 /heap-arrays:0 /traceback /fpp
+/Qdiag-disable:7416 /Qdiag-disable:7025 /Qdiag-disable:5268
+/Fdsrc\libmf6core.a.p\Utilities_compilerversion.F90.pdb
+/module:src\libmf6core.a.p
+/Fosrc/libmf6core.a.p/Utilities_compilerversion.F90.obj /c
+
+System command used to initiate simulation:
+C:\Users\buckl\Documents\Software\mf6.4.1\bin\mf6.exe
+
+MODFLOW was compiled using uniform precision.
+
+Real Variables
+  KIND: 8
+  TINY (smallest non-zero value):    2.225074-308
+  HUGE (largest value):    1.797693+308
+  PRECISION: 15
+  SIZE IN BITS: 64
+
+Integer Variables
+  KIND: 4
+  HUGE (largest value): 2147483647
+  SIZE IN BITS: 32
+
+Long Integer Variables
+  KIND: 8
+  HUGE (largest value): 9223372036854775807
+  SIZE IN BITS: 64
+
+Logical Variables
+  KIND: 4
+  SIZE IN BITS: 32
+ 
+
+ OPENED mfsim.nam
+ FILE TYPE:NAM   UNIT 1001   STATUS:OLD                 
+ FORMAT:FORMATTED              ACCESS:SEQUENTIAL          
+ ACTION:READ                
+
+ # File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+
+ READING SIMULATION OPTIONS
+ END OF SIMULATION OPTIONS
+
+ READING SIMULATION TIMING
+
+ OPENED gwf_bcf2ss.tdis
+ FILE TYPE:TDIS   UNIT 1002   STATUS:OLD                 
+ FORMAT:FORMATTED              ACCESS:SEQUENTIAL          
+ ACTION:READ                
+
+
+ TDIS -- TEMPORAL DISCRETIZATION PACKAGE,
+ VERSION 1 : 11/13/2014 - INPUT READ FROM UNIT 1002
+ # File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+ PROCESSING TDIS OPTIONS
+    SIMULATION TIME UNIT IS DAYS
+ END OF TDIS OPTIONS
+ PROCESSING TDIS DIMENSIONS
+    2 STRESS PERIOD(S) IN SIMULATION
+ END OF TDIS DIMENSIONS
+ PROCESSING TDIS PERIODDATA
+
+
+ STRESS PERIOD     LENGTH       TIME STEPS     MULTIPLIER FOR DELT
+ ----------------------------------------------------------------------------
+        1         1.000000          1                    1.000
+        2         1.000000          1                    1.000
+ END OF TDIS PERIODDATA
+ END OF SIMULATION TIMING
+
+ READING SIMULATION MODELS
+    GWF6 model GWF_BCF2SS will be created as model 1
+ END OF SIMULATION MODELS
+
+ READING SIMULATION EXCHANGES
+ END OF SIMULATION EXCHANGES
+
+ READING SOLUTIONGROUP
+
+ Creating solution: SLN_1           
+
+ OPENED gwf_bcf2ss.ims
+ FILE TYPE:IMS   UNIT 1012   STATUS:OLD                 
+ FORMAT:FORMATTED              ACCESS:SEQUENTIAL          
+ ACTION:READ                
+
+ END OF SIMULATION SOLUTIONGROUP
+
+PROCESSING MODEL CONNECTIONS
+
+ IMS -- ITERATIVE MODEL SOLUTION PACKAGE, VERSION 6, 4/28/2017
+         INPUT READ FROM UNIT 1012
+ # File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+
+ PROCESSING IMS OPTIONS
+ END OF IMS OPTIONS
+
+ PROCESSING IMS NONLINEAR
+ END OF IMS NONLINEAR DATA
+ ***UNDER-RELAXATION WILL NOT BE USED***
+ 
+ ***IMS LINEAR SOLVER WILL BE USED***
+
+ IMSLINEAR -- UNSTRUCTURED LINEAR SOLUTION PACKAGE, VERSION 8, 04/28/2017
+
+ PROCESSING LINEAR DATA
+ END OF LINEAR DATA
+ 
+ An asymmetric matrix will be solved
+
+ OUTER ITERATION CONVERGENCE CRITERION    (DVCLOSE) =    0.100000E-05
+ MAXIMUM NUMBER OF OUTER ITERATIONS        (MXITER) = 500
+ SOLVER PRINTOUT INDEX                     (IPRIMS) = 0
+ NONLINEAR ITERATION METHOD            (NONLINMETH) = 0
+ LINEAR SOLUTION METHOD                   (LINMETH) = 1
+
+       SOLUTION BY THE BICONJUGATE-GRADIENT STABILIZED METHOD
+ ------------------------------------------------------------------
+ MAXIMUM OF 500 CALLS OF SOLUTION ROUTINE
+ MAXIMUM OF 100 INTERNAL ITERATIONS PER CALL TO SOLUTION ROUTINE
+ LINEAR ACCELERATION METHOD            = BCGS      
+ MATRIX PRECONDITIONING TYPE           = MOD. INCOMPLETE LU  
+ MATRIX SCALING APPROACH               = NO SCALING          
+ MATRIX REORDERING APPROACH            = ORIGINAL ORDERING        
+ NUMBER OF ORTHOGONALIZATIONS          = 0
+ HEAD CHANGE CRITERION FOR CLOSURE     =    0.10000E-05
+ RESIDUAL CHANGE CRITERION FOR CLOSURE =    0.10000E-02
+ RESIDUAL CONVERGENCE OPTION           = 1
+ RESIDUAL CONVERGENCE NORM             = INFINITY NORM S 
+ RELAXATION FACTOR                     =    0.97000E+00
+
+
+
+ 
+    Solving:  Stress period:     1    Time step:     1
+ 
+1
+                            STRESS PERIOD NO. 1, LENGTH =   1.000000    
+                            -----------------------------------------------
+                            NUMBER OF TIME STEPS = 1
+                            MULTIPLIER FOR DELT =     1.000
+                            INITIAL TIME STEP SIZE =   1.000000    
+
+
+
+ PSEUDO-TRANSIENT CONTINUATION WILL BE APPLIED TO MODEL 1 ("GWF_BCF2SS") DURING THIS TIME STEP
+ 
+    Solving:  Stress period:     2    Time step:     1
+ 
+1
+                            STRESS PERIOD NO. 2, LENGTH =   1.000000    
+                            -----------------------------------------------
+                            NUMBER OF TIME STEPS = 1
+                            MULTIPLIER FOR DELT =     1.000
+                            INITIAL TIME STEP SIZE =   1.000000    
+
+
+
+ PSEUDO-TRANSIENT CONTINUATION WILL BE APPLIED TO MODEL 1 ("GWF_BCF2SS") DURING THIS TIME STEP
+
+ MEMORY MANAGER TOTAL STORAGE BY DATA TYPE, IN KILOBYTES
+ -------------------------------
+                    ALLOCATED   
+ DATA TYPE           MEMORY     
+ -------------------------------
+ Character        2.4290000     
+ Logical          4.00000000E-02
+ Integer          66.428000     
+ Real             160.47200     
+ -------------------------------
+ Total            229.36900     
+ -------------------------------
+
+
+ 
+ Run end date and time (yyyy/mm/dd hh:mm:ss): 2024/01/13 16:12:58
+ Elapsed run time:  0.042 Seconds
+ Normal termination of simulation.
diff --git a/notebooks/implementing/gwf_bcf2ss/model_p2a/mfsim.nam b/notebooks/implementing/gwf_bcf2ss/model_p2a/mfsim.nam
new file mode 100644
index 0000000..0e67434
--- /dev/null
+++ b/notebooks/implementing/gwf_bcf2ss/model_p2a/mfsim.nam
@@ -0,0 +1,19 @@
+# File generated by Flopy version 3.5.0 on 01/13/2024 at 16:12:15.
+BEGIN options
+END options
+
+BEGIN timing
+  TDIS6  gwf_bcf2ss.tdis
+END timing
+
+BEGIN models
+  gwf6  gwf_bcf2ss.nam  gwf_bcf2ss
+END models
+
+BEGIN exchanges
+END exchanges
+
+BEGIN solutiongroup  1
+  ims6  gwf_bcf2ss.ims  gwf_bcf2ss
+END solutiongroup  1
+
diff --git a/notebooks/learning/gwf_bcf2ss.ipynb b/notebooks/learning/gwf_bcf2ss.ipynb
index 115cd45..3084522 100644
--- a/notebooks/learning/gwf_bcf2ss.ipynb
+++ b/notebooks/learning/gwf_bcf2ss.ipynb
@@ -27,7 +27,10 @@
     "\t\t- 5(b)(iv). River (RIV) Package\n",
     "\t\t- 5(b)(v). Well (WEL) Package\n",
     "\t\t- 5(b)(vi). Recharge (RCH) Package\n",
-    "\t\t- 5(b)(vii). Output Control (OC)"
+    "\t\t- 5(b)(vii). Output Control (OC)\n",
+    "- 6. Breakdown the Newton-Raphson formulation simulation\n",
+    "\t- 6(a). Iterative Model Solution (IMS)\n",
+    "\t- 6(b). Node Property Flow (NPF) Package"
    ]
   },
   {
@@ -2030,6 +2033,169 @@
     "```"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 6. Breakdown the Newton-Raphson formulation simulation\n",
+    "\n",
+    "There are two packages that differ in the Newton-Raphson version: IMS and NPF. The \"NEWTON\" flag is also specified in the model `.nam` file."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 6(a). [Iterative Model Solution (IMS)](https://flopy.readthedocs.io/en/latest/source/flopy.mf6.modflow.mfims.html)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 142,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "package_name = ex-gwf-bcf2ss\n",
+       "filename = ex-gwf-bcf2ss.ims\n",
+       "package_type = ims\n",
+       "model_or_simulation_package = simulation\n",
+       "simulation_name = ex-gwf-bcf2ss\n",
+       "\n",
+       "Block nonlinear\n",
+       "--------------------\n",
+       "outer_dvclose\n",
+       "{internal}\n",
+       "(1e-06)\n",
+       "\n",
+       "outer_maximum\n",
+       "{internal}\n",
+       "(500)\n",
+       "\n",
+       "\n",
+       "Block linear\n",
+       "--------------------\n",
+       "inner_maximum\n",
+       "{internal}\n",
+       "(100)\n",
+       "\n",
+       "inner_dvclose\n",
+       "{internal}\n",
+       "(1e-06)\n",
+       "\n",
+       "rcloserecord\n",
+       "{internal}\n",
+       "(rec.array([('inner_rclose', 0.001, 'strict')],\n",
+       "          dtype=[('inner_rclose_label', 'O'), ('inner_rclose', '<f8'), ('rclose_option', 'O')]))\n",
+       "\n",
+       "linear_acceleration\n",
+       "{internal}\n",
+       "('bicgstab')\n",
+       "\n",
+       "relaxation_factor\n",
+       "{internal}\n",
+       "(0.97)\n",
+       "\n"
+      ]
+     },
+     "execution_count": 142,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "sim_p2a.package_type_dict['ims'][0]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Reverse-engineer class summary!\n",
+    "\n",
+    "```python\n",
+    "outer_dvclose = 1e-6\n",
+    "outer_maximum = 500\n",
+    "inner_maximum = 100\n",
+    "inner_dvclose = 1e-6\n",
+    "rcloserecord = [0.001, \"strict\"]\n",
+    "# preconditioned bi-conjugate gradient stabilized method\n",
+    "linear_acceleration = \"bicgstab\"\n",
+    "relaxation_factor = 0.97\n",
+    "filename = \"ex-gwf-bcf2ss.ims\"\n",
+    "pname = \"ex-gwf-bcf2ss\"\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 6(b). [Node Property Flow (NPF) Package](https://flopy.readthedocs.io/en/latest/source/flopy.mf6.modflow.mfgwfnpf.html)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 143,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "package_name = npf\n",
+       "filename = ex-gwf-bcf2ss.npf\n",
+       "package_type = npf\n",
+       "model_or_simulation_package = model\n",
+       "model_name = ex-gwf-bcf2ss\n",
+       "\n",
+       "Block options\n",
+       "--------------------\n",
+       "save_specific_discharge\n",
+       "{internal}\n",
+       "(True)\n",
+       "\n",
+       "\n",
+       "Block griddata\n",
+       "--------------------\n",
+       "icelltype\n",
+       "Layer_1{constant 1}\n",
+       "Layer_2{constant 0}\n",
+       "\n",
+       "k\n",
+       "Layer_1{constant 10.0}\n",
+       "Layer_2{constant 5.0}\n",
+       "\n",
+       "k33\n",
+       "{constant 0.1}\n",
+       "\n"
+      ]
+     },
+     "execution_count": 143,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "model_p2a.package_dict['npf']"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Reverse-engineer class summary!\n",
+    "\n",
+    "```python\n",
+    "save_specific_discharge = True\n",
+    "icelltype = [1, 0]\n",
+    "k = [10.0, 5.0]\n",
+    "k33 = 0.1\n",
+    "filename = \"ex-gwf-bcf2ss.npf\"\n",
+    "pname = \"npf\"\n",
+    "```"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,