Workshop_7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4491806c-4cae-43b5-b5cf-cf15d04d94ff",
   "metadata": {},
   "source": [
    "# MODPATH and ZoneBudget\n",
    "Hi and welcome to Workshop 7. In this session we will look at two commonly used tools ZoneBudget and MODPATH that are often applied in conjunction with MODFLOW type models. The former is used for assessing fluxes from all sinks and sources in a section of a model designated by a zone. The latter is used for particle tracing, which can be a more appropriate alternative to conservative solute transport for certain modelling questions. We will demonstrate basic application of each with an applied use case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08cfde0c-339a-4ee3-93f1-437fdd5ff96b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import shutil\n",
    "import platform\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import flopy\n",
    "from flopy.discretization import VertexGrid\n",
    "from flopy.utils import Raster\n",
    "from flopy.utils import GridIntersect\n",
    "from flopy.utils.gridgen import Gridgen\n",
    "gridgen_exe = \"gridgen\"\n",
    "if platform.system() in \"Windows\":\n",
    "    gridgen_exe += \".exe\"\n",
    "gridgen_exe = flopy.which(\"gridgen\")\n",
    "if gridgen_exe is None:\n",
    "    msg = (\n",
    "        \"Warning, gridgen is not in your path. \"\n",
    "        \"When you create the griden object you will need to \"\n",
    "        \"provide a full path to the gridgen binary executable.\"\n",
    "    )\n",
    "    print(msg)\n",
    "else:\n",
    "    print(\"gridgen executable was found at: {}\".format(gridgen_exe))\n",
    "\n",
    "print(f\"Pandas version = {pd.__version__}\")\n",
    "print(f\"Numpy version = {np.__version__}\")\n",
    "print(f\"Flopy version = {flopy.__version__}\")\n",
    "print(f\"Matplotlib version = {mpl.__version__}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "863668f2-a0e3-4970-b87d-1ff6cbfc2994",
   "metadata": {},
   "outputs": [],
   "source": [
    "ws7 = os.path.join('workshop_7') # here we are making a path not creating the folder\n",
    "gis_f = os.path.join(ws7,'GIS') # creating a sub-directory path for our gis input/output\n",
    "model_f = os.path.join(ws7,'model') # creating a sub-directory path for our model input/output\n",
    "plots_f = os.path.join(ws7,'plots') # creating a sub-directory path for our plots\n",
    "for path in [ws7,gis_f,model_f,plots_f]:\n",
    "    if os.path.exists(path): # here we are asking if the path exists on the computer. \n",
    "        shutil.rmtree(path)# if it does exist, delete it and all the files in it\n",
    "        os.mkdir(path) # then remake it\n",
    "    else:\n",
    "        os.mkdir(path) # if it doesn't exist then make the folder"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3e4a7cd-d3bd-40ba-b076-4b0e745efa0b",
   "metadata": {},
   "source": [
    "# Import a model buid from scratch\n",
    "This takes some time. This is the same model that we've used previously in many other tutorials. It includes the MAW wells that we configured in the previous session."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "166f57f8-3ccb-43bc-bbcb-6f0fa7937f4e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from helpers import ws7_model1\n",
    "sim, gwf = ws7_model1(ws7,gis_f,model_f,plots_f)\n",
    "# ignore the warnings"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a70fa0b-bcfc-4ed6-ae54-3d18555221a0",
   "metadata": {},
   "source": [
    "# If the model already exists just load an rerun\n",
    "Just a reminder that if you have a pre-existing model you can load it using something similar to the code block below. Make sure that you get the model name correct or it will throw an error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26cab773-5419-4cae-ad78-2c9f6fc5277c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sim = flopy.mf6.MFSimulation.load(sim_name=\"MySim\",sim_ws=model_f)\n",
    "# gwf = sim.get_model(model_name='flow')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89b2a6bb-75eb-4218-bcdf-59b02d1ce865",
   "metadata": {},
   "source": [
    "# Running the imported simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32dfd97c-e7fd-4430-a1fd-37b3028ab451",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d68ba12-3fe8-46ad-8673-75c35a8eabed",
   "metadata": {},
   "source": [
    "# Saving SS heads for use later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3987e536-8c8b-4b4f-823e-c5ca4ae7381f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save initial heads for later use\n",
    "headfile = os.path.join(model_f,\"flow.hds\")\n",
    "hds = flopy.utils.binaryfile.HeadFile(headfile)\n",
    "h = hds.get_data((0,0))\n",
    "np.savetxt(\"iheads_array.txt\",h[0]) # note this is not saved to the model folder because if I choose to rerun the script from scrtach it will be deleted."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4fd9330-34de-4cd2-afe2-ed7ba502abba",
   "metadata": {},
   "source": [
    "# Adjusting tdis\n",
    "Here we take advantage of the record array features to quickly make adjustments to the tdis input file via a Flopy object. All I want to do is increase the number of time steps in each stress period. Note the object needs to be created from the imported simulation. Our import returned both the sim object and the model object so we have these available to work from although we could have just returned the sim object alone in this case. First we create the tdis object then check the period data. Then we access the data directly via the .array attribute (it is a record array after all)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ed6ccb8-4627-453d-a69b-803ac76ee6fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "tdis = sim.get_package(\"tdis\")\n",
    "tdis.perioddata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63c851d5-d8d2-458d-9224-ad5bafd4ff30",
   "metadata": {},
   "outputs": [],
   "source": [
    "tdis.perioddata.array.nstp[1:]=4\n",
    "tdis.perioddata"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9be6ef2-abcf-4921-b679-0ab6644d1226",
   "metadata": {},
   "source": [
    "# MODPATH and start_date_time\n",
    "Modpath does not like the option start_date_time in the tdis object so we have to remove this as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "072fcc0a-45ef-4583-8497-4b0f57cf9229",
   "metadata": {},
   "outputs": [],
   "source": [
    "tdis.start_date_time=None\n",
    "tdis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba5e5e95-64f7-44e3-ad9d-980f3ef988f0",
   "metadata": {},
   "source": [
    "# Change OC to output \"ALL\"\n",
    "It is preferable to output a higher resolution of time intervals when using particle tracking especially when doing transient tracing through layers with different porosity. The manual for MODPATH goes into a bit more detail as to why this is the case. In the code block below there are several different dictionaries for different output setups. Take note of which ones are being used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74d9ae9b-757f-4788-b7e9-5643c855f560",
   "metadata": {},
   "outputs": [],
   "source": [
    "test = list(range(1,sim.tdis.nper.data)) # this range represents the monthly stress periods before recovery\n",
    "hs_keys = [0,*test] # SSkey = 0, then all monthly SP\n",
    "\n",
    "# use this as the saverecord data in the oc object if we don't want a budget file\n",
    "h_rec = {key:[(\"HEAD\",\"LAST\")] for key in hs_keys}\n",
    "\n",
    "# use this as the saverecord data for combined head plus budget for zonebudget run\n",
    "zbud_rec = {key:[(\"BUDGET\",\"LAST\"),(\"HEAD\",\"LAST\")] for key in hs_keys}\n",
    "\n",
    "# use this as the saverecord data for combined head plus budget for MODPATH run\n",
    "mpbud_rec = {key:[(\"BUDGET\",\"ALL\"),(\"HEAD\",\"ALL\")] for key in hs_keys}\n",
    "\n",
    "#use this as the printrecord data for budget printing to list file\n",
    "b_rec = {key:[(\"BUDGET\",\"LAST\")] for key in hs_keys}\n",
    "\n",
    "oc = flopy.mf6.ModflowGwfoc(\n",
    "    gwf,\n",
    "    pname=\"oc\",\n",
    "    budget_filerecord=\"{}.cbb\".format(gwf.name),\n",
    "    head_filerecord=\"{}.hds\".format(gwf.name),\n",
    "    headprintrecord=[(\"COLUMNS\", 10, \"WIDTH\", 15, \"DIGITS\", 6, \"GENERAL\")],\n",
    "    saverecord=mpbud_rec,\n",
    "    printrecord=b_rec,\n",
    "    budgetcsv_filerecord=\"myflowbudget.csv\"\n",
    ")\n",
    "oc.write()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1348a7aa-7ae0-4228-9126-520840af271b",
   "metadata": {},
   "source": [
    "# Change initial conditions\n",
    "Here we are just setting it to use the heads we saved earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d94f6c4-c1fa-48f5-8946-90adab362a5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# using hdata from before\n",
    "mg=gwf.modelgrid\n",
    "ihd_array=np.ones_like(mg.botm)\n",
    "last_ssheads = np.loadtxt(\"iheads_array.txt\")\n",
    "ihd_array[:]=last_ssheads\n",
    "\n",
    "ic = flopy.mf6.ModflowGwfic(\n",
    "    gwf, pname=\"ic\", strt=ihd_array, filename=\"flow.ic\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f43945b0-3073-4844-8d0d-e087e50cdb27",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.write_simulation()\n",
    "sim.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecfbaa33-5095-4850-806d-57a8f171dbb8",
   "metadata": {},
   "source": [
    "# Setting up particles to trace wells\n",
    "We are going to setup particles for back-tracing to be released from multiple locations around model cells that the MAW package uses. We don't want all the particles to release from the same location in each cell but rather different positions around the model cell. First step is to get the node numbers. To be compatible with both DISV and DISU MODPATH uses global nodes, that is, it treats DISV nodes the same as DISU. Recall that DISV cells have id (layer, nodenumber) while DISU only have node number. The global node value for a DISV cell in layer two is given by the node  number above it in layer 1 plus the number of cells in each layer. For layer 3 this would be the node number above it in layer 1 plus 2 times the number of cells in each layer. So the first step is to get the well nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a088ce4-5a83-441e-90c9-09f7f931a9b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get well nodes\n",
    "maw = gwf.get_package('maw')\n",
    "well_cells = maw.connectiondata.array.cellid\n",
    "well_cells\n",
    "# we only want the lowest node that the wells are connected to\n",
    "l2_well_cells = [x for x in well_cells if x[0]==1]\n",
    "l1_only_wells = [x for x in well_cells if not (1,x[1]) in l2_well_cells]\n",
    "lowest = [*l1_only_wells,*l2_well_cells]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d43ce83-5aea-4131-82b6-fe0223e41000",
   "metadata": {},
   "outputs": [],
   "source": [
    "lowest"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b02e45d-2fc7-4044-b0a2-e0e749fe46e4",
   "metadata": {},
   "source": [
    "# DISV and NCPL\n",
    "NCPL is a property of a DISV type grid and is an integer value of the number of cells in a layer. We can use this to get the global node values for nodes in other layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "136942bc-f2ad-4f7d-8f67-21935b897efe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the global nodes using \n",
    "nodews = [gwf.disv.ncpl.array * x + y for (x,y) in lowest]\n",
    "# global nodes\n",
    "nodews"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dc3534b-89f5-42ef-b44d-fefe61650149",
   "metadata": {},
   "source": [
    "# Setting up positions around a cell for particle release\n",
    "We don't want all cells starting at the exact same location in a cell so we will offset them from the centre of the cell. If you consider the node being central to the x y and z axis of a given cell then it's position within the cell is half the distance of the x, y and z dimensions of the cell. The block below sets x and y positions for particles around a given cell. There are 16 positions, four per model cell edge but all at the node elevation (z = 0.5)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b38cd5b3-2d3b-47ca-8db6-1f88c0fdfdea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use this to get localx, localy and localz around each node\n",
    "import numpy as np\n",
    "ls = [round(x, 3) for x in [0.125,0.375,0.625,0.875]]\n",
    "ls2 = [[x,y,0.5] for x in [0.0,1.0] for y in ls]\n",
    "ls3 = [[y,x,0.5] for x in [0.0,1.0] for y in ls]\n",
    "ls4 = [*ls2,*ls3]\n",
    "pcoord = np.asarray(ls4)\n",
    "pcoord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9743c76d-28ef-4eb0-8f79-eed32f5ab732",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pcoord.shape, pcoord.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1302b090-3806-4fe1-b260-771ffe5fd7c5",
   "metadata": {},
   "source": [
    "# Pairing up particels with nodes\n",
    "So we will be releasing 16 particles from each model cell given by the cells local x', local y and local z positions described in the array we created above. Now we need to pair each of those model cell local positions with the actual node numbers. This means we will need 16 repeating node numbers, one to pair with each of the model cell local positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13fb884c-ba60-412c-bf6a-b163964091bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "plocs = [nodew for nodew in nodews for i in range(pcoord.shape[0])]\n",
    "plocs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac711642-108f-43e3-9c9e-9b141b97e0c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(plocs))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95cc41d6-d125-4762-a984-1e999678c15f",
   "metadata": {},
   "source": [
    "# Repeating the local positions\n",
    "For each model cell that we are going to release particles from there needs to be a record of the local position that the particle will start at. This means that the array we created with the 16 positions will need to be repeated the same number of times for each model cell. What we are aiming for is that the total number of particles and the total number of local x,y,z values are the same. In this case 240 particles will be released from 240 different locations. Although some of these locations have a common node number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "297a6253-777d-47c6-aebd-71e1337e73f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "pcoords_all = pcoord.copy()\n",
    "for i in range(len(nodews)-1):\n",
    "    pcoords_all = np.concatenate((pcoords_all,pcoord)) # useing concatenate to basically repeat the array\n",
    "print(pcoords_all.shape)\n",
    "pcoords_all"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a7cd741-32dc-4949-8c1b-da2f36757fa6",
   "metadata": {},
   "source": [
    "# The Particle data\n",
    "So all the effort thus far has been to obtain the input for a modpath particle data object. It requires particle locations provided by the nodenumbers (plocs) and the local node or model cell offsets (localx, localy, localz)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dd4be26-df14-4549-8e2a-91c66623eba5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create particle data\n",
    "pa = flopy.modpath.ParticleData(\n",
    "    plocs,\n",
    "    structured=False, # True when using l,r,c False otherwise\n",
    "    localx=pcoords_all[:, 0],\n",
    "    localy=pcoords_all[:, 1],\n",
    "    localz=pcoords_all[:, 2],\n",
    "    drape=0, # has to do with wet and dry cell behaviour\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dace8fc-888a-4a4f-b6ed-397253f4f227",
   "metadata": {},
   "source": [
    "# Create a particle group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d1d8b57-7b9b-4a6b-b223-9c14cd920026",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create backward particle group\n",
    "fpth = f\"{gwf.name}_mp1.sloc\" # using the model name\n",
    "\n",
    "pga = flopy.modpath.ParticleGroup(\n",
    "    particlegroupname=\"backtrace1\", # can be anything you want\n",
    "    particledata=pa, # pass in the particle data object\n",
    "    filename=fpth # use the filename we created above\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5279cad-4026-4325-a1d4-9ed989622b85",
   "metadata": {},
   "source": [
    "# Create a MODPATH model object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "097d5d6d-1de5-4df7-84ee-b48f4f3bbcf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create modpath files\n",
    "\n",
    "# create a modpath model name\n",
    "mpmod_nam = f\"{gwf.name}_mp1\"\n",
    "\n",
    "# create the modpath model object\n",
    "mp = flopy.modpath.Modpath7(\n",
    "    modelname=mpmod_nam, \n",
    "    flowmodel=gwf, \n",
    "    exe_name=\"mp7\", \n",
    "    model_ws=model_f\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91f9da96-8486-45e0-a6c0-8408b1af709c",
   "metadata": {},
   "source": [
    "# Configure a porosity array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fa552c2-1771-49cf-9493-9a34553c1092",
   "metadata": {},
   "outputs": [],
   "source": [
    "# build a simple porosity array\n",
    "por = np.ones_like(mg.botm)\n",
    "por[0] = 0.2\n",
    "por[1] = 0.2\n",
    "por[2] = 0.05\n",
    "# define the porosity note\n",
    "mpbas = flopy.modpath.Modpath7Bas(mp, porosity=por)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d010efef-fbfa-42be-be86-ec7c5ad0ad93",
   "metadata": {},
   "source": [
    "# Create a MODPATH simualtion object\n",
    "The simulation level object is where the majority of the simulation behaviour information is entered. These need careful consideration because it can have a big impact on your traces and time series information. The concept of reference time is important to understand because this has flow on effects for other settings as well. Knowing the difference between a weak and strong sink or source is also important as is the effect that porosity has on transient trace models. Our simulation does not run for very long and we are only going to consider the transient part of the simulation. Nevertheless it does serve as a useful demonstration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c2239a0-4877-40d9-b74b-0d09ae923273",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the simulation level modpath object\n",
    "flopy.modpath.Modpath7Sim(\n",
    "    mp,\n",
    "    simulationtype=\"combined\",\n",
    "    trackingdirection=\"backward\",\n",
    "    weaksinkoption=\"stop_at\", # stopping at weak sinks\n",
    "    weaksourceoption=\"stop_at\", # stopping at weak sources\n",
    "    referencetime=367.0, # this is the end of the model runtime because we are doing backtracing\n",
    "    stoptimeoption=\"specified\",\n",
    "    stoptime= 366.0, # before the start of the steady-state stress period\n",
    "    timepointdata=[52, 7.0], # release particles every 7 days for a total of 365 days \n",
    "    particlegroups=pga,\n",
    ")\n",
    "\n",
    "# write modpath datasets\n",
    "mp.write_input()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6205324-da0a-4de3-990c-a6fd995e8b3b",
   "metadata": {},
   "source": [
    "# Run it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dacb8cda-c49d-4c50-8280-ca0e45348da5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run modpath\n",
    "success, buff = mp.run_model(silent=True, report=True)\n",
    "assert success, \"mp7 failed to run\"\n",
    "for line in buff:\n",
    "    print(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02361b83-fc94-4bb1-abe9-65ea1e6cff13",
   "metadata": {},
   "source": [
    "# Create a pathline object\n",
    "Here we are pulling information from one of the output files and creating an object that will have methods which we can use for further processing. This contains the particle trace information. The time series information associated with particle locations is in a different output file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8846e129-1b2a-44de-9206-c847daad3ccf",
   "metadata": {},
   "outputs": [],
   "source": [
    "fpth = os.path.join(model_f, f\"{mpmod_nam}.mppth\")\n",
    "p = flopy.utils.PathlineFile(fpth)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3433573-b77b-4a63-a0a6-6b4784613ed8",
   "metadata": {},
   "source": [
    "# Get the data from the pathline object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28aa897e-9782-4a14-ae35-a104dc45434f",
   "metadata": {},
   "outputs": [],
   "source": [
    "p0=p.get_alldata()\n",
    "len(p0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "216921b5-843d-44fd-bc91-51c18edd98f3",
   "metadata": {},
   "source": [
    "# Create a particle timeseries object\n",
    "Here we pull the information from the timeseries output file to create the timeseries object and we get all the data similar to what we did previously for the traces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c030ec70-0cbd-4cda-834f-209b7439ca23",
   "metadata": {},
   "outputs": [],
   "source": [
    "fpth = os.path.join(model_f,f\"{mpmod_nam}.timeseries\")\n",
    "ts = flopy.utils.TimeseriesFile(fpth)\n",
    "ts0 = ts.get_alldata()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52de505b-21ce-458f-87de-b96607977dc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "wc_l1_nodes = [x[1] for x in well_cells]\n",
    "pl1 = p.get_destination_pathline_data(wc_l1_nodes, to_recarray=True)\n",
    "tl1 = ts.get_destination_timeseries_data(wc_l1_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a5caa2d-ef8e-410a-b83f-c0cc0d548150",
   "metadata": {},
   "outputs": [],
   "source": [
    "disv = gwf.get_package(\"dis\")\n",
    "nlay = disv.nlay.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08d66df7-e941-49a6-a9ca-026fc8455d6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "test = np.zeros_like((mg.botm), dtype=int)\n",
    "test = test.flatten()\n",
    "test[nodews] = 1\n",
    "test = test.reshape((mg.botm.shape))\n",
    "ibd = np.ma.masked_equal(test, 0)\n",
    "colors = [\"green\", \"orange\", \"red\"]\n",
    "cmap = mpl.colors.ListedColormap([\"r\",\"g\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a615fb42-f078-4809-a821-a68fa2ad20a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, axs = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 3, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(24, 8)) # setting figure size Note had to increase the x dimension here\n",
    "\n",
    "for i,ax in enumerate(axs):\n",
    "    pmv = flopy.plot.PlotMapView(modelgrid=mg, layer = i, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "    pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "    v = pmv.plot_array(ibd, cmap=cmap, edgecolor=\"gray\", ax=ax)\n",
    "    pmv.plot_pathline(p0, layer=i, colors=\"blue\", lw=0.75, alpha=0.3)\n",
    "    pmv.plot_timeseries(ts0, layer=i, marker=\"o\", lw=0, color=colors[i], alpha=0.1)\n",
    "    ax.set_title(f'Model Layer {i+1}')\n",
    "    ax.set_xlabel('Eastings')\n",
    "    ax.set_ylabel('Northings')\n",
    "    ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis\n",
    "plt.tight_layout() # reduces the whitespace around the figure so it takes up less space in a document"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cf2ca1d-a404-481c-96c4-99b02c0705a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, axs = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 3, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(12, 8),\n",
    "                        constrained_layout=True) # setting figure size Note had to increase the x dimension here\n",
    "\n",
    "for i,ax in enumerate(axs):\n",
    "    pmv = flopy.plot.PlotMapView(modelgrid=mg, layer = i, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "    pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "    v = pmv.plot_array(ibd, cmap=cmap, edgecolor=\"gray\", ax=ax)\n",
    "    pmv.plot_pathline(p0, layer=i, colors=\"blue\", lw=0.75, alpha=0.3)\n",
    "    pmv.plot_timeseries(ts0, layer=i, marker=\"o\", lw=0, color=colors[i], alpha=0.1)\n",
    "    ax.set_title(f'Model Layer {i+1}')\n",
    "    ax.set_xlabel('Eastings')\n",
    "    ax.set_ylabel('Northings')\n",
    "    ax.set_xlim(755750,757500)\n",
    "    ax.set_ylim(963500,967000)\n",
    "    ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "448a7ec8-7cca-4456-9b3c-3f13d2b565b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, axs = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 3, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(12, 8)) # setting figure size Note had to increase the x dimension here\n",
    "\n",
    "for i,ax in enumerate(axs):\n",
    "    pmv = flopy.plot.PlotMapView(modelgrid=mg, layer = i, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "    pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "    v = pmv.plot_array(ibd, cmap=cmap, edgecolor=\"gray\", ax=ax)\n",
    "    pmv.plot_pathline(p0, layer=i, colors=\"blue\", lw=0.75, alpha=0.3)\n",
    "    pmv.plot_timeseries(ts0, layer=i, marker=\"o\", lw=0, color=colors[i], alpha=0.1)\n",
    "    ax.set_title(f'Model Layer {i+1}')\n",
    "    ax.set_xlim(756500,757000)\n",
    "    ax.set_ylim(965000,965500)\n",
    "    ax.set_xlabel('Eastings')\n",
    "    ax.set_ylabel('Northings')\n",
    "    ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis\n",
    "plt.tight_layout() # reduces the whitespace around the figure so it takes up less space in a document"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ffec36b-a662-40e9-9dce-015a023348a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, axs = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 3, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(12, 8)) # setting figure size Note had to increase the x dimension here\n",
    "\n",
    "for i,ax in enumerate(axs):\n",
    "    pmv = flopy.plot.PlotMapView(modelgrid=mg, layer = i, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "    pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "    v = pmv.plot_array(ibd, cmap=cmap, edgecolor=\"gray\", ax=ax)\n",
    "    pmv.plot_pathline(p0, layer=i, colors=\"blue\", lw=0.75, alpha=0.3)\n",
    "    #pmv.plot_timeseries(ts0, layer=i, marker=\"o\", lw=0, color=colors[i], alpha=0.1)\n",
    "    ax.set_title(f'Model Layer {i+1}')\n",
    "    ax.set_xlim(756500,757000)\n",
    "    ax.set_ylim(965000,965500)\n",
    "    ax.set_xlabel('Eastings')\n",
    "    ax.set_ylabel('Northings')\n",
    "    ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis\n",
    "plt.tight_layout() # reduces the whitespace around the figure so it takes up less space in a document"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a63cbeac-2753-4cb7-be7f-017e32fb2a92",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hmm not much happening there"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ef4a271-53e4-4bc7-91c9-b64c622a0c38",
   "metadata": {},
   "outputs": [],
   "source": [
    "sp_dict = {}\n",
    "sp_dict[0] = [[x,mg.zcellcenters[x]+5] for x in lowest]\n",
    "chd = flopy.mf6.ModflowGwfchd(gwf,stress_period_data=sp_dict)\n",
    "gwf.remove_package(\"maw\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bd9ebc9-c3d7-4310-872b-439516359c8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "gwf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50149dbc-da1e-4e27-a95f-816316bc8fbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "tdis.perioddata.array.perlen[0]=36500.0\n",
    "tdis.perioddata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf90ca33-24a8-4d8e-87a3-7ee283f4b68b",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.write_simulation()\n",
    "sim.run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ef7cfd2-fcc6-42d6-8956-09244e2b734e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the simulation level modpath object\n",
    "flopy.modpath.Modpath7Sim(\n",
    "    mp,\n",
    "    simulationtype=\"combined\",\n",
    "    trackingdirection=\"backward\",\n",
    "    weaksinkoption=\"pass_through\", # stopping at weak sinks\n",
    "    weaksourceoption=\"pass_through\", # stopping at weak sources\n",
    "    referencetime=1.0, # this is the end of the model runtime because we are doing backtracing\n",
    "    stoptimeoption=\"extend\",\n",
    "    timepointdata=[100,365.0], \n",
    "    particlegroups=pga,\n",
    ")\n",
    "\n",
    "# write modpath datasets\n",
    "mp.write_input()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf0a8a1e-a1ab-4e97-a448-343caefe9d28",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run modpath\n",
    "success, buff = mp.run_model(silent=True, report=True)\n",
    "assert success, \"mp7 failed to run\"\n",
    "for line in buff:\n",
    "    print(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05eb9d70-e16f-4a9e-95a5-646d10913e71",
   "metadata": {},
   "outputs": [],
   "source": [
    "fpth = os.path.join(model_f, f\"{mpmod_nam}.mppth\")\n",
    "p = flopy.utils.PathlineFile(fpth)\n",
    "p0=p.get_alldata()\n",
    "\n",
    "fpth = os.path.join(model_f,f\"{mpmod_nam}.timeseries\")\n",
    "ts = flopy.utils.TimeseriesFile(fpth)\n",
    "ts0 = ts.get_alldata()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e78e6c04-5392-42b6-b32d-7eba09f8413c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, axs = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 3, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(24, 8)) # setting figure size Note had to increase the x dimension here\n",
    "\n",
    "for i,ax in enumerate(axs):\n",
    "    pmv = flopy.plot.PlotMapView(modelgrid=mg, layer = i, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "    pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "    v = pmv.plot_array(ibd, cmap=cmap, edgecolor=\"gray\", ax=ax)\n",
    "    pmv.plot_pathline(p0, layer=i, colors=\"blue\", lw=0.75, alpha=0.3)\n",
    "    pmv.plot_timeseries(ts0, layer=i, marker=\"o\", lw=0, color=colors[i], alpha=0.1)\n",
    "    ax.set_title(f'Model Layer {i+1}')\n",
    "    ax.set_xlabel('Eastings')\n",
    "    ax.set_ylabel('Northings')\n",
    "    ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis\n",
    "plt.tight_layout() # reduces the whitespace around the figure so it takes up less space in a document"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6062474-d01e-42a9-86b0-e82c5f062773",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, axs = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 3, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(12, 8),\n",
    "                        constrained_layout=True) # setting figure size Note had to increase the x dimension here\n",
    "\n",
    "for i,ax in enumerate(axs):\n",
    "    pmv = flopy.plot.PlotMapView(modelgrid=mg, layer = i, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "    pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "    v = pmv.plot_array(ibd, cmap=cmap, edgecolor=\"gray\", ax=ax)\n",
    "    pmv.plot_pathline(p0, layer=i, colors=\"blue\", lw=0.75, alpha=0.3)\n",
    "    #pmv.plot_timeseries(ts0, layer=i, marker=\"o\", lw=0, color=colors[i], alpha=0.1)\n",
    "    ax.set_title(f'Model Layer {i+1}')\n",
    "    ax.set_xlabel('Eastings')\n",
    "    ax.set_ylabel('Northings')\n",
    "    ax.set_xlim(755750,757500)\n",
    "    ax.set_ylim(963500,967000)\n",
    "    ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f8052e2-e72a-4167-a875-c9fe06fb7f2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, axs = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 3, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(12, 8)) # setting figure size Note had to increase the x dimension here\n",
    "\n",
    "for i,ax in enumerate(axs):\n",
    "    pmv = flopy.plot.PlotMapView(modelgrid=mg, layer = i, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "    pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "    v = pmv.plot_array(ibd, cmap=cmap, edgecolor=\"gray\", ax=ax)\n",
    "    pmv.plot_pathline(p0, layer=i, colors=\"blue\", lw=0.75, alpha=0.3)\n",
    "    #pmv.plot_timeseries(ts0, layer=i, marker=\"o\", lw=0, color=colors[i], alpha=0.1)\n",
    "    ax.set_title(f'Model Layer {i+1}')\n",
    "    ax.set_xlim(756500,757000)\n",
    "    ax.set_ylim(965000,965500)\n",
    "    ax.set_xlabel('Eastings')\n",
    "    ax.set_ylabel('Northings')\n",
    "    ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis\n",
    "plt.tight_layout() # reduces the whitespace around the figure so it takes up less space in a document"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21816983-0a6f-44d7-ab75-61b7d143882c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new modpath model name\n",
    "mpmod_nam2 = f\"{gwf.name}_mp2\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccf808a6-b06d-4d53-873a-b3ccff619bca",
   "metadata": {},
   "outputs": [],
   "source": [
    "facedata = flopy.modpath.FaceDataType(\n",
    "    drape=0,\n",
    "    verticaldivisions1=10,\n",
    "    horizontaldivisions1=10,\n",
    "    verticaldivisions2=10,\n",
    "    horizontaldivisions2=10,\n",
    "    verticaldivisions3=10,\n",
    "    horizontaldivisions3=10,\n",
    "    verticaldivisions4=10,\n",
    "    horizontaldivisions4=10,\n",
    "    rowdivisions5=0,\n",
    "    columndivisions5=0,\n",
    "    rowdivisions6=4,\n",
    "    columndivisions6=4,\n",
    ")\n",
    "\n",
    "pb = flopy.modpath.NodeParticleData(subdivisiondata=facedata, nodes=nodews)\n",
    "# create forward particle group\n",
    "fpth = f\"{mpmod_nam2}.sloc\"\n",
    "pgb = flopy.modpath.ParticleGroupNodeTemplate(\n",
    "    particlegroupname=\"backtrace2\", particledata=pb, filename=fpth\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb16d85b-1b61-4753-98bd-9981b5cdeacd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create modpath files\n",
    "mp = flopy.modpath.Modpath7(\n",
    "    modelname=mpmod_nam2, \n",
    "    flowmodel=gwf, \n",
    "    exe_name=\"mp7\", \n",
    "    model_ws=model_f\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b5a58f6-8d5d-45ea-a4b8-830469948ee4",
   "metadata": {},
   "outputs": [],
   "source": [
    "flopy.modpath.Modpath7Bas(mp, porosity=por)\n",
    "\n",
    "# generally always use this for endpoint analysis\n",
    "flopy.modpath.Modpath7Sim(\n",
    "    mp,\n",
    "    simulationtype=\"endpoint\",\n",
    "    trackingdirection=\"backward\",\n",
    "    weaksinkoption=\"pass_through\",\n",
    "    weaksourceoption=\"pass_through\",\n",
    "    referencetime=0.0,\n",
    "    stoptimeoption=\"extend\",\n",
    "    particlegroups=pgb,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "634f1cd5-948b-4225-bdd3-772357edecae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write modpath datasets\n",
    "mp.write_input()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "895f1292-8867-44af-811f-e13346981b20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run modpath\n",
    "success, buff = mp.run_model(silent=True, report=True)\n",
    "assert success, \"mp7 failed to run\"\n",
    "for line in buff:\n",
    "    print(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f4211c3-8abe-4aef-9a6a-7d3ddbe71f40",
   "metadata": {},
   "outputs": [],
   "source": [
    "fpth = os.path.join(model_f, f\"{mpmod_nam2}.mpend\")\n",
    "e = flopy.utils.EndpointFile(fpth)\n",
    "e0 = e.get_alldata()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49531e21-a198-40fd-9473-7bebce74281f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, ax = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 1, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(8, 8), # setting figure size Note had to increase the x dimension here\n",
    "                        constrained_layout=True)\n",
    "                        \n",
    "\n",
    "pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "v = pmv.plot_array(ibd, cmap=cmap, edgecolor=\"gray\", ax=ax)\n",
    "pmv.plot_endpoint(e0, direction=\"ending\", colorbar=True, shrink=0.5)\n",
    "ax.set_title(f'Model Layer 1 Endpoint Analysis')\n",
    "ax.set_xlabel('Eastings')\n",
    "ax.set_ylabel('Northings')\n",
    "ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2213134b-deb3-4626-914f-e7381edbebdb",
   "metadata": {},
   "source": [
    "# Zonebudget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a103d13-b0b4-45b1-81d4-022118392f2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make anothe rmodel folder\n",
    "model2_f = os.path.join(ws7,'model2') # creating a sub-directory path for our model input/output\n",
    "if os.path.exists(model2_f): # here we are asking if the path exists on the computer. \n",
    "    shutil.rmtree(model2_f)# if it does exist, delete it and all the files in it\n",
    "    os.mkdir(model2_f) # then remake it\n",
    "else:\n",
    "    os.mkdir(model2_f) # if it doesn't exist then make the folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5d92461-2b99-4e87-9031-91b294f53326",
   "metadata": {},
   "outputs": [],
   "source": [
    "m2files_path = os.path.join(\"files\",\"zbudmod\") \n",
    "flist = [x for x in os.listdir(m2files_path)]\n",
    "for file in flist:\n",
    "    shutil.copyfile(os.path.join(m2files_path,file),os.path.join(model2_f,file))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd7fdb88-5a4f-46b9-a615-6dafd90afcfe",
   "metadata": {},
   "source": [
    "# Load the new simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "231894b1-18cc-4325-bec1-5e201df66cf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = flopy.mf6.MFSimulation.load(sim_name=\"MySim2\",sim_ws=model2_f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9120b48c-8b30-43d4-9efa-9d7b2b638d92",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80f11ab9-0dde-4c70-8daf-7a21bff14f38",
   "metadata": {},
   "source": [
    "# Get the model object from a loaded model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b71cc877-dd97-4687-a6c4-9e45758bcbeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "gwf = sim.get_model(model_name=\"flow\")\n",
    "mg=gwf.modelgrid\n",
    "gwf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68473f67-3b85-4975-8fc2-97f490ac1721",
   "metadata": {},
   "source": [
    "# Get some boundary conditions from laoded model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0442380-9de4-4e82-9457-e0671ca76bf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the perimeter bounds\n",
    "ghb_ds = gwf.get_package(name=\"ghb_ds\")\n",
    "ghb_dv = gwf.get_package(name=\"ghb_dv\")\n",
    "ghb_dd = gwf.get_package(name=\"ghb_dd\")\n",
    "ghb_ea = gwf.get_package(name=\"ghb_ea\")\n",
    "ghb_sw = gwf.get_package(name=\"ghb_sw\")\n",
    "ghb_so = gwf.get_package(name=\"ghb_so\")\n",
    "ghb_ls = [ghb_ds,ghb_dv,ghb_dd,ghb_ea,ghb_sw,ghb_so]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "186316c5-1817-47d1-9740-5bfe2ea0f83d",
   "metadata": {},
   "source": [
    "# Take a look at model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "640a1445-2bd3-4db0-aec6-4ab1bbb52f63",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, ax = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 1, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(7, 5), # setting figure size Note had to increase the x dimension here\n",
    "                        constrained_layout=True)\n",
    "                        \n",
    "\n",
    "pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "for bc in ghb_ls:\n",
    "    pmv.plot_bc(package=bc, color='red')\n",
    "# Plot a shapefile of a polygon representing the water control district\n",
    "shp = os.path.join(model2_f, \"wd.shp\") # path to the shapefile of a line I created in GIS\n",
    "# plot the cross section line\n",
    "patch_collection = pmv.plot_shapefile(shp, \n",
    "                                      radius=0, \n",
    "                                      lw=3,\n",
    "                                      alpha = 0.1,\n",
    "                                      edgecolor=\"blue\", \n",
    "                                      facecolor=\"blue\")\n",
    "ax.set_title(f'Model2 Layer 1')\n",
    "ax.set_xlabel('Eastings')\n",
    "ax.set_ylabel('Northings')\n",
    "ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb62c4ec-4dd9-4185-8c2c-8e27aa0e7a85",
   "metadata": {},
   "source": [
    "# Extract heads for processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f3bc012-8908-4150-b31e-3d75f6c9e7f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the heads file\n",
    "headfile = os.path.join(model2_f,\"flow.hds\")\n",
    "hds = flopy.utils.binaryfile.HeadFile(headfile)\n",
    "# get the heads from the steady-state stress period\n",
    "h0 = hds.get_data((0,0))\n",
    "end_pumping = (11,23)\n",
    "# get the heads from the final stress period- Note there is only one step in each stress period.\n",
    "h = hds.get_data(end_pumping)\n",
    "# get our contour levels using linspace\n",
    "levels = np.linspace(np.min(h),np.max(h),20)\n",
    "print(levels)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77cf9812-cde7-4650-8764-fde8da9e6381",
   "metadata": {},
   "source": [
    "# Plot SS heads and contours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "875babd9-6105-4be5-9bbf-1909dd5071c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, ax = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 1, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(7, 5), # setting figure size Note had to increase the x dimension here\n",
    "                        constrained_layout=True)\n",
    "                        \n",
    "\n",
    "pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "pc = pmv.plot_array(h0)\n",
    "colorbar = plt.colorbar(pc, aspect=30, shrink= 0.8) # create a colorbar\n",
    "contour_set = pmv.contour_array(h0, levels=levels, colors='w',linewidths = 0.75) # plot some contours using our levels\n",
    "ax.clabel(contour_set, fmt='%.1f', colors='w', fontsize=8) # change the appearance of the contours\n",
    "ax.set_title(f'Model2 Layer 1 initial SS heads (mAHD)')\n",
    "ax.set_xlabel('Eastings')\n",
    "ax.set_ylabel('Northings')\n",
    "ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9589b8af-a858-45e0-a13a-54467baec399",
   "metadata": {},
   "source": [
    "# Process heads for drawdown"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14047fb0-3206-467f-8e09-4bc91dc3d0f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the drawdown\n",
    "ddn = h0-h\n",
    "ddn[ddn>50.0]=50.0\n",
    "ddn[ddn<0.1]=0.0\n",
    "levels = [0.1, 1.0, 5.0, 10.0, 20.0, 50.0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54294d96-e90b-4fde-831d-7d9b6f62038e",
   "metadata": {},
   "source": [
    "# Plot our drawdown"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcfad69e-2624-4a35-a9ff-42a685b19912",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, ax = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 1, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(7, 5), # setting figure size Note had to increase the x dimension here\n",
    "                        constrained_layout=True)\n",
    "                        \n",
    "\n",
    "pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "pc = pmv.plot_array(ddn,masked_values=[0.0]) # note that I mask where there is no drawdown\n",
    "colorbar = plt.colorbar(pc, aspect=30, shrink= 0.8) # create a colorbar\n",
    "contour_set = pmv.contour_array(ddn, levels=levels, colors='w',linewidths = 0.75) # plot some contours using our levels\n",
    "ax.clabel(contour_set, fmt='%.1f', colors='w', fontsize=8) # change the appearance of the contours\n",
    "ax.set_title(f'Model2 Layer 1 end of pumping drawdown (m)')\n",
    "ax.set_xlabel('Eastings')\n",
    "ax.set_ylabel('Northings')\n",
    "ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d0eee1e-deeb-4c98-981a-1962a2122573",
   "metadata": {},
   "source": [
    "# Get a zone array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2663594-b414-43ed-96fe-8a8a922708cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "def get_bnodes(shpfyl): # works with single and multiple line strings\n",
    "    ix = GridIntersect(mg, method=\"vertex\") # build our intersection object we will be applyitng this to a vertex grid only\n",
    "    poly = gpd.read_file(shpfyl).geometry # read in the shapefile\n",
    "    if len(poly)==1: # if the feature only has one entry\n",
    "        return(ix.intersect(poly[0]).cellids) # note this is an array\n",
    "    else: \n",
    "        ls = [] # if the feature has mutiple items i.e. points or multi-line strings\n",
    "        for item in poly: # loop through the different geometries of each item\n",
    "            nums = ix.intersect(item).cellids\n",
    "            ls = [*ls,*nums]\n",
    "        return(np.asarray(ls)) # # this will also be an array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d8b48fc-6549-4684-b838-00aa755f618a",
   "metadata": {},
   "outputs": [],
   "source": [
    "zdf = get_bnodes(shp)\n",
    "zdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e52fe398-8efb-461b-a383-32c685c05880",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we need to make a zone array\n",
    "zones = np.ones_like(mg.botm)*2\n",
    "zones[0][zdf.astype('int')]=1.0\n",
    "zones"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "117b9d04-ee6c-4cdd-a512-b0cc7b7eb1cb",
   "metadata": {},
   "source": [
    "# Plot our Zones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d75741ea-cce5-4378-be82-ac3f8ac3cf2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using Matplotlib subplots to create a mutiple axes figure\n",
    "fig, ax = plt.subplots(nrows = 1, # I want one row\n",
    "                        ncols = 1, # I want three columns\n",
    "                        sharey=True, # they can share the Y-axis to save space\n",
    "                        subplot_kw={'aspect':'equal'}, # Each sub_plot must have aspect = equal to keep x and y ratios consistent\n",
    "                        figsize=(7, 5), # setting figure size Note had to increase the x dimension here\n",
    "                        constrained_layout=True)\n",
    "                        \n",
    "\n",
    "pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) #creating a mapview object assigning a specific layer and specifying the axes\n",
    "pmv.plot_grid(ax=ax, lw=0.3, color=\"black\") # plot the grid on the axes using linewidths = 0.3 and black lines\n",
    "pc = pmv.plot_array(zones) # note that I mask where there is no drawdown\n",
    "colorbar = plt.colorbar(pc, aspect=30, shrink= 0.8) # create a colorbar\n",
    "ax.set_title(f'Model2 Layer 1 Zones')\n",
    "ax.set_xlabel('Eastings')\n",
    "ax.set_ylabel('Northings')\n",
    "ax.ticklabel_format(style='plain') #  gets rid of the exponent offsets on the axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "276504e9-652c-423d-850f-6d9c69f25283",
   "metadata": {},
   "outputs": [],
   "source": [
    "from flopy.utils import ZoneBudget\n",
    "zb6_exe = \"zbud6.exe\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96db6ca8-4929-418d-b764-cce7a7d157fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "aliases = {1: \"WDWCD\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d22c9f5-9e90-4489-b98d-6bf82234fb5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "zonbud = gwf.output.zonebudget(zones.astype(\"int\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "410ec99f-b578-4220-ae65-55b5f39c5288",
   "metadata": {},
   "outputs": [],
   "source": [
    "grbnam = gwf.name+'.disv.grb'\n",
    "zonbud.add_package(\"grb\",grbnam)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54d99906-804e-4e89-89eb-98ad97e08b09",
   "metadata": {},
   "source": [
    "# Write input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2190909-acd6-441f-a0d5-45bdee4068a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "zonbud.write_input()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b8028a1-0abe-496c-8861-5a1ff7a647f5",
   "metadata": {},
   "source": [
    "# Run the zonebudget model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5273412-76a1-4465-bcd5-b467e74343b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "zonbud.run_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fea8205-0402-4111-9027-9d823eb2e9b0",
   "metadata": {},
   "source": [
    "# Get the budget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f0aef09-4149-4b49-a650-44b180f75a21",
   "metadata": {},
   "outputs": [],
   "source": [
    "zonbud.get_budget()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1792acdb-da53-4477-8690-92a0b09487ed",
   "metadata": {},
   "source": [
    "# Complete model budget to and from zones\n",
    "Make sure you understand how to read this layout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c08c3fda-e0f7-44d9-ab57-889dcee8831a",
   "metadata": {},
   "outputs": [],
   "source": [
    "zdf = zonbud.get_dataframes(net=True, pivot=True)\n",
    "zdf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47c43c2e-c24f-404c-aac6-07a2feb1c1ff",
   "metadata": {},
   "source": [
    "# Get the y-axis values as flux from zone 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49c89ec4-cb62-44b0-8267-d618e5ec896c",
   "metadata": {},
   "outputs": [],
   "source": [
    "wd_flux = zonbud.get_dataframes(names='FROM_ZONE_1', pivot=True)\n",
    "final_df = wd_flux[wd_flux['zone']==2]\n",
    "final_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fde3655-13ea-4c67-b002-8e8a4f7dce76",
   "metadata": {},
   "source": [
    "# Get x axis values as years"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31ad7868-ba6c-4e3e-9873-5ae1ca5c0775",
   "metadata": {},
   "outputs": [],
   "source": [
    "times = final_df[\"totim\"]\n",
    "times = np.round(times/365.25,2)\n",
    "times"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccb923f6-0f5d-45de-bb27-9c97aec6ca62",
   "metadata": {},
   "source": [
    "# Plot the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81cb5eee-2d49-4e2e-8e1e-a6aaf43a3125",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(times,final_df['FROM_ZONE_1']/1000)\n",
    "ax.set_ylabel(\"Flux ($ML/d$)\")\n",
    "ax.set_xlabel(\"Project Years\")\n",
    "ax.set_title(\"Flux from controlled district \")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29492367-cd3d-467f-a71b-c12fbefd8280",
   "metadata": {},
   "source": [
    "# Thats all Folks!\n",
    "Hopefully we have some time for Q/A and any requests for something to cover in the final session. If not then we will have a look at some other handy features Flopy makes available to you for modelling purposes."
   ]
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}