Workshop_3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "72c198ae-4f61-4492-95c3-479a4c2aac0a",
   "metadata": {},
   "source": [
    "# Workshop 3\r\n",
    "Hi and welcome to the third workshop. This session will focus on boundary conditions, specifically list type boundary conditions. Some boundary conditions also allow for specifying arrays but we will cover that in the next session along with other array style input packages. You will also be introduced to use and configuration of time series files plus use and configuration of observation files. Additionally, we will introduce specifying auxiliary variables as either multipliers or for transport modelling.  \r\n",
    "\r\n",
    "As before please make sure you have the MF6io.pdf document open while working through this notebook because it will show you the model input file structure for the specific boundary/package that you are creating and provides a more complete reference of the different options availabl.\r\n",
    "."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3804a7d-3e8a-429b-a985-80869403a3c8",
   "metadata": {},
   "source": [
    "# Imports first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9bbe270-7652-4eb0-850b-2dc9ac5cab30",
   "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\n",
    "import matplotlib.pyplot as plt\n",
    "import flopy\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 = {matplotlib.__version__}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "956b521d-3850-44d0-b8ab-991340c75776",
   "metadata": {},
   "source": [
    "# Build a simple steady-state structured model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "934d9e56-90c7-48b5-b79b-0a6073f4ffee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup some folders\n",
    "ws3 = os.path.join('workshop_3') # here we are making a path not creating the folder\n",
    "gis_f = os.path.join(ws3,'GIS') # creating a sub-directory path for our gis input/output\n",
    "model_f = os.path.join(ws3,'model') # creating a sub-directory path for our model input/output\n",
    "plots_f = os.path.join(ws3,'plots') # creating a sub-directory path for our plots\n",
    "for path in [ws3,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\n",
    "\n",
    "#Build some basic model objects\n",
    "simulation_name = 'MySim' # setting a name for our simulation\n",
    "sim = flopy.mf6.MFSimulation(sim_name=simulation_name, sim_ws=model_f, exe_name='mf6') # for a no frills mfsim.nam file\n",
    "tdis = flopy.mf6.ModflowTdis(sim) # for the MySim.tdis file, recall this sets a single steady-state stress period\n",
    "ims = flopy.mf6.ModflowIms(sim) # for the MySim.ims file, recall this sets up IMS as 'SIMPLE'\n",
    "model_name = 'flow' # setting a name for our model\n",
    "gwf = flopy.mf6.ModflowGwf(sim, modelname=model_name, save_flows=True) # for the flow.nam file\n",
    "\n",
    "# Recall that there are default values applied to whatever you don't specify. \n",
    "# What esle should have been specified? Check your mf6io.pdf document.\n",
    "# Make sure you are aware of the default entries that Flopy uses when you don't specify something for a package.\n",
    "# REMEBER\n",
    "# *** as the modeller it is your responsibility to check your input files ***\n",
    "nrow = 20\n",
    "ncol = 20\n",
    "dis = flopy.mf6.ModflowGwfdis(gwf, nrow=nrow, ncol=ncol) # why not just specify the values in here directly? Why use ncol=ncol?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f15445e-5263-40bd-b544-07f5443f092a",
   "metadata": {},
   "source": [
    "# Taking a look at things"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31bb087d-2f7b-4114-9a66-2bfa07de5ae9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets get a modelgrid object and have a look at the grid\n",
    "mg = gwf.modelgrid\n",
    "\n",
    "# create a figure to save to our plots folder\n",
    "fig,ax = plt.subplots(figsize=(5,5)) # we are creating a figure object here so that we can dictate size\n",
    "mg.plot(ax=ax)\n",
    "ax.ticklabel_format(style='plain') # test what happens if you don't use this switch\n",
    "ax.set_title('DIS Model Grid')\n",
    "ax.set_xlabel('Eastings')\n",
    "ax.set_ylabel('Northings')\n",
    "# lets save it to our plots folder\n",
    "figname = os.path.join(plots_f,'DIS_model_grid.png') # note use of the plots folder path.\n",
    "# If you want to change the file format then change the extension from .png to pdf or just do both\n",
    "fig.savefig(figname,dpi=300)\n",
    "figname = os.path.join(plots_f,'DIS_model_grid.pdf')\n",
    "fig.savefig(figname,dpi=300) \n",
    "\n",
    "# While we are here we can export it as a shapefile too.\n",
    "dis.export(os.path.join(gis_f,'dis.shp'))\n",
    "# ignore the warnings"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afcbe7d6-f174-4ade-b791-de3d3dddef02",
   "metadata": {},
   "source": [
    "# The stress period data dictionary\n",
    "This is arguably the most important part of setting up boundary conditions. Once you grasp how to use and build the information needed in the stress period dictionary, the setup of a basic boundary conditions becomes reasonably straight forward. Whatever variables you see in the MF6io.pdf document between BEGIN PERIOD X and END PERIOD X will have to included in the dictionary. A dictionary is not always necessary as will be demonstrated but learning how to use them will be of benefit when it comes to setting up more complex versions of the same boundary. Boundary conditions are specific to a model so the model object must exist and be provided. Check the options block of whatever boundary condition you want to setup in the MF6io.pdf. Most of the settings in the options block are set with simple entries like option = \"True\", option= None, option = string, option = list_of_strings or option = integer etc. Check the Flopy defaults for each option and make sure you understand what it means to turn different settings on or off."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa319d5c-64c6-457a-8b32-1750e5f19db2",
   "metadata": {},
   "source": [
    "# Active or Inactive\n",
    "You need to pay special attention to how the specific boundary you are configuring is made active or inactive. This varies for some boundaries and is noted in the MF6io.pdf document, so please be sure to read the specifics for the boundary you are working on. For the boundaries we will demonstrate in this workshop they are active from the first stress period you designate and will remain active until explicitly deactivated in a subsequent stress period. If you don't designate a stress period in which a boundary condition should first activate then Flopy assumes you want the boundary active throughout the entire simulation. You should always check your package input file before running your model anyway to ensure that everything is correct."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e72cfb9f-0e42-40a3-a9d4-0cc9585c081b",
   "metadata": {},
   "source": [
    "# Assign a constant head\n",
    "To begin with we will assign a first type boundary because it is relatively simple and then we will demonstrate all the different boundary condition features offered in MF6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f99d18ba-47e5-4124-8e4b-2dbeee56f071",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assigning a constant head bounadry to specific model cells to be active all the time.\n",
    "\n",
    "chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=[[(0, 0, 0), 10.],\n",
    "                                                       [(0, 19, 19), 8.]])\n",
    "\n",
    "# which model cells were assigned? what value were they given?\n",
    "# recall we have a structured model so what do the numbers in the tuples represent?\n",
    "# If this were a DISV grid, how would this input change?\n",
    "# What if it was DISU?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a57a81a1-568b-42a5-b51f-b7ef60c44c83",
   "metadata": {},
   "source": [
    "# Writing your boundary condition package\n",
    "If you want to be able to write your input file on its own then you need to name the object. What I am referring to here is the \"chd = \" part. You can instantiate the object without doing this and it will still work. We have been assigning names to the objects all along because it allows us to easily access the built in methods for that object. You've seen how we write a simulation before in Workshop 1. Writing the simulation will write all the model files you've created so far but here we only want the CHD package file we just created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "715f3163-aad3-4197-8e06-f46dd755045c",
   "metadata": {},
   "outputs": [],
   "source": [
    "chd.write()\n",
    "# note that it will determine what value needs to be assigned to the maxbound variable \n",
    "# based on the information you provide in your stress_period_data.\n",
    "\n",
    "# Look for the file in your model folder and open it using your referred text editor or use the box below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33e356b7-f321-4c7f-ae1f-6179da2f536d",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,f\"{model_name}.chd\"))]\n",
    "# note the name of the input file will be the model name followed by the package Ftype acronym. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67b47665-b547-4fd5-9794-4eb6ab305764",
   "metadata": {},
   "source": [
    "# Plotting the boundary to check it\n",
    "It is always a good idea to visually check that the cells you intended to set the boundary on are correct. You can do this by plotting the model grid and also plotting a specific boundary condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a55326fe-f365-4646-92ce-82c141abbdb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the same method we applied previously in Workshop 2 to plot the modelgrid using a MapView object\n",
    "# remeber the sequence of steps\n",
    "# 1. Get modelgrid\n",
    "# 2. Create figure instance\n",
    "# 3. Create Mapview object\n",
    "# 4. Tell the Mapview object what you want to plot\n",
    "# 5. Modify your figure and axes to your preferred settings\n",
    "# Note the line used to plot the bounadry condition\n",
    "\n",
    "# Lets get a modelgrid object and have a look at the grid\n",
    "mg = gwf.modelgrid\n",
    "\n",
    "# create a figure instance to save to our plots folder\n",
    "fig,ax = plt.subplots(figsize=(7,5)) # we are creating a figure object here so that we can dictate size note there are mutiple ways to do this\n",
    "\n",
    "# Create the mapview object of the modelgrid and tell it to use the axes we created with our figure\n",
    "pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) # note by not specifying a layer here it will assume layer 1\n",
    "\n",
    "# say that you want the Mapview object to plot the grid with linewidth = 0.3 and colour the lines black\n",
    "pmv.plot_grid(ax=ax, lw=0.3, color=\"black\")\n",
    "\n",
    "# Now I also want to plot the boundary condition so I specify the package object to plot\n",
    "pmv.plot_bc(package=chd)\n",
    "\n",
    "\n",
    "ax.set_title('CHD Boundary Cells')\n",
    "ax.set_xlabel('Eastings')\n",
    "ax.set_ylabel('Northings')\n",
    "# lets save it to our plots folder\n",
    "figname = os.path.join(plots_f,'CHD_bounadry_cells.png') # note use of the plots folder path.\n",
    "# If you want to change the file format then change the extension from .png to pdf or just do both\n",
    "fig.savefig(figname,dpi=300)\n",
    "figname = os.path.join(plots_f,'CHD_bounadry_cells.pdf')\n",
    "fig.savefig(figname,dpi=300) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1c0ddf3-1119-4e57-9ac5-30d379df18ee",
   "metadata": {},
   "source": [
    "There are other options available when plotting a boundary condition that we will explore later. This includes specifying the layer you wish to view if you have a multi-layer model, specifying a stress period if you have boundary conditions that vary between active and inactive during your simulation, changing the colour of the plotted boundary and if you want to plot all boundaries."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44aaa757-0199-46f5-944e-a0f8a5da3648",
   "metadata": {},
   "source": [
    "# Selecting cells for a boundary condition\r\n",
    "The first example demonstrates how to select individual cells in a structured grid but we may also want to assign a boundary condition along the edge of a model grid. Selecting model cells to assign a boundary condition to is half the battle so having different methods for selecting cells is going o help you a lot. We'll stick with our structured grid for now and look at selecting cells using some list comprehensions. We are going to assign our boundaries down the left and right side of our model as opposed to the corners. This means we need to build our stress period dictionary a little bit differently\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "394e080f-ca56-4138-80cd-61efdeb1466f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets first get rid of our existing chd object\n",
    "del chd \n",
    "\n",
    "# We only have one layer so if using l,r,c it means that we always have l=1 (or 0 for zero based)\n",
    "l=0\n",
    "# If we are on the left boundary then c is always 1 (or 0 for zero based) \n",
    "c = 0\n",
    "# this means that the only variable we need to loop over is r\n",
    "left_bound_nodes = [(l,r,c) for r in range(nrow)]\n",
    "# lets take a look at it\n",
    "left_bound_nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "821df2db-3d3a-42c4-9199-893bf4a3ea40",
   "metadata": {},
   "source": [
    "So that sorts out the nodes but the stress period dictionary needs more information than that. We also need to specify the value of the constant head in those cells. Each entry in the period block of the input file needs to provided to Flopy as a list with the node id first then the rest of the boundary information, which is the constant head value in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a8fee51-7701-46ae-a757-26efa5e02759",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we just add on the rest of the information that we want to include in the list comprehension and make each loop provide a list.\n",
    "left_bound_stress_period_block = [[(l,r,c),10.] for r in range(nrow)]\n",
    "left_bound_stress_period_block\n",
    "# note what you produced is a list of lists evidenced by the leading and trailing double brackets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ca4dcda-6fad-4f86-8f17-f335b4a01894",
   "metadata": {},
   "outputs": [],
   "source": [
    "# but we also need it for the right boundary so lets do the same thing here\n",
    "# c on the right boundary is the number of columns minus one for zero base\n",
    "c = ncol-1\n",
    "right_bound_stress_period_block = [[(l,r,c),8.] for r in range(nrow)]\n",
    "right_bound_stress_period_block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e8224ff-646a-4927-9e16-2fee6be03a00",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We need our left and right boundaries incorporated into a single list for the stress period dictionary so lets unpack them\n",
    "chd_period_list = [*right_bound_stress_period_block,*left_bound_stress_period_block]\n",
    "chd_period_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8981d76-e501-4157-af51-e0fe414aeec2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can build our new CHD bounadry but this time we provide our list as an object name\n",
    "chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_period_list)\n",
    "# write it\n",
    "chd.write()\n",
    "# Look at the file\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,f\"{model_name}.chd\"))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c7c9814-4d3f-49da-8916-b52c5d15fecc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Then plot it to have a better look\n",
    "\n",
    "# Lets get a modelgrid object and have a look at the grid\n",
    "mg = gwf.modelgrid\n",
    "\n",
    "# create a figure instance to save to our plots folder\n",
    "fig,ax = plt.subplots(figsize=(7,5)) # we are creating a figure object here so that we can dictate size note there are mutiple ways to do this\n",
    "\n",
    "# Create the mapview object of the modelgrid and tell it to use the axes we created with our figure\n",
    "pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) # note by not specifying a layer here it will assume layer 1\n",
    "\n",
    "# say that you want the Mapview object to plot the grid with linewidth = 0.3 and colour the lines black\n",
    "pmv.plot_grid(ax=ax, lw=0.3, color=\"black\")\n",
    "\n",
    "# Now I also want to plot the boundary condition so I specify the package object to plot\n",
    "pmv.plot_bc(package=chd)\n",
    "\n",
    "\n",
    "ax.set_title('CHD Boundary Cells')\n",
    "ax.set_xlabel('Eastings')\n",
    "ax.set_ylabel('Northings')\n",
    "# lets save it to our plots folder\n",
    "figname = os.path.join(plots_f,'CHD_bounadry_cells_LR.png') # note use of the plots folder path.\n",
    "# If you want to change the file format then change the extension from .png to pdf or just do both\n",
    "fig.savefig(figname,dpi=300)\n",
    "figname = os.path.join(plots_f,'CHD_bounadry_cells_LR.pdf')\n",
    "fig.savefig(figname,dpi=300) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17b0ae2e-32bd-403b-9809-1b2a920687c9",
   "metadata": {},
   "source": [
    "# But wait there is a problem here\r\n",
    "Take a look at the model object. There are two boundary packages assigned to the model object even though we deleted the first chd object and used the same object name. MF6 allows multiples of the same boundary type to be used. Previous versions of model only allowed one instance of each boundary package. Look at the default naming convention used by Flopy when you have multiple instances of the same package. You can overrule this through the options block, which we will do shortly but first we need to fix our extra chd package option.\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd681fde-b9c8-42d4-8a85-8dbec52e463a",
   "metadata": {},
   "outputs": [],
   "source": [
    "gwf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fff30e84-0d13-4bb0-8e1f-86eb56350682",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(gwf.get_package_list())\n",
    "gwf.remove_package(\"CHD_0\")\n",
    "print(gwf.get_package_list())\n",
    "# but which one did we keep?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8da0735e-06b7-4b90-b03c-82662e367fe8",
   "metadata": {},
   "source": [
    "But which one did we keep? Check by plotting to be sure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ef7ed56-ae5d-462d-99d3-c7ca98595f59",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Lets get a modelgrid object and have a look at the grid\n",
    "mg = gwf.modelgrid\n",
    "\n",
    "# create a figure instance to save to our plots folder\n",
    "fig,ax = plt.subplots(figsize=(7,5)) # we are creating a figure object here so that we can dictate size note there are mutiple ways to do this\n",
    "\n",
    "# Create the mapview object of the modelgrid and tell it to use the axes we created with our figure\n",
    "pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) # note by not specifying a layer here it will assume layer 1\n",
    "\n",
    "# say that you want the Mapview object to plot the grid with linewidth = 0.3 and colour the lines black\n",
    "pmv.plot_grid(ax=ax, lw=0.3, color=\"black\")\n",
    "\n",
    "# Now I also want to plot the boundary condition so I specify the package object to plot\n",
    "pmv.plot_bc(package=chd)\n",
    "\n",
    "\n",
    "ax.set_title('CHD Boundary Cells')\n",
    "ax.set_xlabel('Eastings')\n",
    "ax.set_ylabel('Northings')\n",
    "# lets save it to our plots folder\n",
    "figname = os.path.join(plots_f,'CHD_bounadry_cells_LR.png') # note use of the plots folder path.\n",
    "# If you want to change the file format then change the extension from .png to pdf or just do both\n",
    "fig.savefig(figname,dpi=300)\n",
    "figname = os.path.join(plots_f,'CHD_bounadry_cells_LR.pdf')\n",
    "fig.savefig(figname,dpi=300) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03d2a0ab-1284-446d-b7fc-786404d4c00c",
   "metadata": {},
   "source": [
    "Lets delete the old input files and re-write our package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b83267c-cddb-487e-a3a3-d26221e136a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = [os.remove(os.path.join(model_f,x)) for x in os.listdir(model_f)]\n",
    "# for x in os.listdir(model_f) means give me all the names of the files in the model folder\n",
    "# os.remove(os.path.join(model_f,x)) means remove the file at the path given by os.path.join(model_f,x) \n",
    "\n",
    "print(os.listdir(model_f)) # should be an empty list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3d23157-de12-4d43-a40b-64f4e23d15cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now write our model file\n",
    "chd.write()\n",
    "\n",
    "# but which file name will it use now?\n",
    "print(os.listdir(model_f))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df68911e-7e01-4992-8845-d1dd3a554eb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at the file\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,'flow_0.chd'))]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f957e3c-8127-4b83-8576-7d562e7f034e",
   "metadata": {},
   "source": [
    "# Run a model\r\n",
    "The following import will provide a function for you to use so that you can run a model using what you have setup so far. This is similar to the introductory model on the Flopy GitHub site. It will also do some post processing for you so that we can view results quickly and focus on the boundary conditions.\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a56b9c2d-35d3-488e-94e6-6cb227c4ab8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from helpers import ws3_mod1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "167e8a89-7be6-4d52-82a5-b23ca55ab147",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clear out all existing modle files\n",
    "_ = [os.remove(os.path.join(model_f,x)) for x in os.listdir(model_f)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de666118-c155-4e42-9c4d-798d46cfe310",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the function by passing in your sim and gwf objects\n",
    "ws3_mod1(sim,gwf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "686b1b17-cb30-417b-bb40-8a29412faf7f",
   "metadata": {},
   "source": [
    "Let's try adding a well rouglhy at the centre of our model that extracts 1L/s (864 m3/d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48f3f3a2-34aa-4e93-978d-9350e611fc59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Take a look at the MF6io.pdf document to see what is needed in the well stress period data block\n",
    "# then follow the same procedure we did previously\n",
    "wel1 = [(0,9,9), -10.0]\n",
    "wel_period_list = [wel1] # remember that the stress_period data is a list of lists\n",
    "wel = flopy.mf6.ModflowGwfwel(gwf, stress_period_data=wel_period_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75fb459b-5a63-48ea-a9fd-50fdeff6da43",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check our model object to make sure the well package is there\n",
    "gwf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d3ed6f1-39df-43de-85d1-8524f49ad9ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clear out all existing modle files\n",
    "_ = [os.remove(os.path.join(model_f,x)) for x in os.listdir(model_f)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29cf24dd-71a2-4385-9937-f64413cc95e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run our model again\n",
    "ws3_mod1(sim,gwf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "384de507-24ed-4900-b24e-ab803ac4d5a9",
   "metadata": {},
   "source": [
    "# What if I want to know the fluxes associated with the boundaries?\r\n",
    "You can get this from the list file but this workshop is on boundary conditions so we'll get straight from there. First, we will have to assign our boundaries some names and then setup an observation package to go with each of them. Now we will start to use some of the options in the options block. Again take a look at the MF6io.pdf document to check what you can do with each boundary condition. We will start by removing our CHD and WEL packages.\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27b04630-e36f-46c6-bd51-1155258d3bbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(gwf.get_package_list())\n",
    "gwf.remove_package(\"CHD\") \n",
    "gwf.remove_package(\"WEL\")\n",
    "print(gwf.get_package_list())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e305a9fe-61f2-463e-8783-cc8372f00981",
   "metadata": {},
   "source": [
    "# Why don't I have to remove the other packages as well?\r\n",
    "The remaining packages will be overwritten because you can only have one of them assigned to a model. You can have several boundary conditions of the same type in each model so you need to delete the ones that you don't want or you may accidentally double up. Again it is your responsibility to check your input files before running your model.\r\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a799d878-adc3-45b7-98b8-ec8cb8fd8562",
   "metadata": {},
   "source": [
    "# Adding boundary names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a023398a-320c-48f0-875f-6dbf3d014eaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now our stress period data list will have to include some more information so we need to rebuild it\n",
    "c = 0\n",
    "# this is what we had before [[(l,r,c),10.]\n",
    "left_bound_stress_period_block = [[(l,r,c),10.0,\"left_chd\"] for r in range(nrow)] # All we did was add in a string for this bounadry name\n",
    "# now the right one\n",
    "c = ncol-1\n",
    "right_bound_stress_period_block = [[(l,r,c),8.0,\"right_chd\"] for r in range(nrow)]\n",
    "chd_period_list = [*left_bound_stress_period_block,*right_bound_stress_period_block]\n",
    "chd_period_list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6508f865-c02a-4cd6-a966-9672b26518ee",
   "metadata": {},
   "source": [
    "# Setting up a boundary observation package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8b6e8ce-87cf-4c7d-8dc3-b77818d77dae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# build obs data dictionary\n",
    "chd_obs = {\n",
    "    \"chd_flux.csv\": [                    # This is the user specified name of the file that the boundary fluxes will be stored in after the model has run\n",
    "        (\"Flux_left\", \"CHD\", \"left_chd\"), # A tuple with what will be the column heading in the output csv file, the code for the type of observation , then the boundary name\n",
    "        (\"Flux_right\", \"CHD\", \"right_chd\"), # A tuple with what will be the column heading in the output csv file, the code for the type of observation , then the boundary name\n",
    "    ],\n",
    "}\n",
    "# So what you have is a dictionary with the output filename as the key\n",
    "# note how the boundary name is used to efectivley group a bunch of boundary cells together\n",
    "# We could have broken the boundary up into more sections with more names if wanted\n",
    "# The code informing the observation package what it needs to observe differs for each package\n",
    "# We will look through these in the workshop\n",
    "chd_obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38e6ff7e-acf7-46e6-be3a-5dc349842fee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clear out all existing modle files\n",
    "_ = [os.remove(os.path.join(model_f,x)) for x in os.listdir(model_f)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1b91a4a-5a6f-4556-bcfd-22326778c5a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "chd = flopy.mf6.ModflowGwfchd(gwf, \n",
    "                              boundnames=True, # Tell it we will be using boundary names\n",
    "                              observations=chd_obs, # Tell it we will be using an observation package\n",
    "                              stress_period_data=chd_period_list)\n",
    "# write the boundary and its observation package\n",
    "chd.write()\n",
    "chd.obs.write()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54f50e11-ad5c-411c-beed-9bddf108c2f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,'flow.chd'))]\n",
    "# Note the options have been added in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae09eaa6-7ee7-4742-9c4c-2a2b5b211d93",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,'flow.chd.obs'))]\n",
    "# this is what the observation file looks like\n",
    "# check that it makes sense according to the mf6io.pdf document"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46c69831-0d28-4b30-9e41-f98fccfa117b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we will do it for the well but in a slightly different way\n",
    "wel1 = [(0,9,9), -10.0,\"Mywell\"] # adding in a boundary name\n",
    "wel_period_list = [wel1] # remember that the stress_period data is a list of lists\n",
    "wel = flopy.mf6.ModflowGwfwel(gwf, boundnames=True, stress_period_data=wel_period_list) \n",
    "# added the bound and built the object but didn't setup an observation package\n",
    "# Lets do it now using an initialization method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e49bb6cc-414d-40d3-ac19-bb3dc90b33ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup the observation dictionary the same as before\n",
    "wel_obs = {\n",
    "    \"wel_flux.csv\": [\n",
    "        (\"Well rate\", \"WEL\", \"Mywell\")]\n",
    "}\n",
    "wel.obs.initialize(continuous=wel_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8c63b4e-e50a-410a-b76d-d00c1bfd193b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write the boundary and its observation package\n",
    "wel.write()\n",
    "wel.obs.write()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c77afe6-0cce-4ae5-821f-dacd84ba8ccc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets run our model to see what happens\n",
    "ws3_mod1(sim,gwf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6519b09-2179-4e6a-a9cf-ecaaa7e7489b",
   "metadata": {},
   "source": [
    "# Lets try a transient model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e7e0376-be11-4f43-8d02-404e2a5df6bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we'll remove the existing tdis package from sim object\n",
    "sim.remove_package(\"TDIS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1432d391-48f1-4e18-a20f-8591a8eade88",
   "metadata": {},
   "outputs": [],
   "source": [
    "del tdis # delete the existing tdis object\n",
    "# now we build a new one\n",
    "tdis = flopy.mf6.ModflowTdis(sim,\n",
    "                             time_units='days',\n",
    "                             start_date_time=\"2023-12-31T00:00:00\",\n",
    "                             nper=2,\n",
    "                             perioddata=[(1,1,1),(365.25,52,1)]) # steady-state initial stress period\n",
    "tdis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc8f54f7-0a61-45f1-bfa1-1cf5ac082056",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import a transient version of the same model we used before\n",
    "from helpers import ws3_mod1_trans"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00c27bb2-7620-4037-97f5-c3da618a9f8f",
   "metadata": {},
   "source": [
    "# Changing to a transient stress period data dictionary\r\n",
    "So far we have been using a list of lists to create our stress period data and this is fine for boundary conditions that don't vary throughout your simulation. But for this transient simulation we don't want our well active in the steady-state stress period. We want it to start in the second stress period.\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "181e3fe2-84ba-4f6d-8621-23bbe3af3a1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove the existing wel package from our simulation\n",
    "gwf.remove_package(\"WEL\")\n",
    "# Now we build a new one starting with the stress period data DICTIONARY\n",
    "# We still use the list of lists but now it is stored inside a dictionary \n",
    "# with the key representing the stress period number\n",
    "\n",
    "wel_period_dict = {} # start an empty dictionary\n",
    "wel1 = [(0,9,9), -10.0,\"Mywell\"] # adding in a boundary name\n",
    "wel_period_list = [wel1] # remember that the stress_period data is a list of lists\n",
    "wel_period_dict[1] = wel_period_list # note the dictioanry key is 1, this is stress period two (zero base)\n",
    "\n",
    "wel = flopy.mf6.ModflowGwfwel(gwf, boundnames=True, stress_period_data=wel_period_dict) # Note we changed this to the dictionary \n",
    "# added the bound and built the object but didn't setup an observation package\n",
    "# Lets do it now using an initialization method\n",
    "\n",
    "# setup the observation dictionary the same as before\n",
    "wel_obs = {\n",
    "    \"wel_flux.csv\": [\n",
    "        (\"Well rate\", \"WEL\", \"Mywell\")]\n",
    "}\n",
    "wel.obs.initialize(continuous=wel_obs)\n",
    "# check the model object to make sure we haven't doubled up on boundaries\n",
    "gwf.get_package_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "def0b68a-95f8-4b26-853f-a238f9134a49",
   "metadata": {},
   "outputs": [],
   "source": [
    "ws3_mod1_trans(sim,gwf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bf89486-69fd-46df-97c8-8bad4da4a29a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets look at the well file\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,'flow.wel'))]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cedfbf92-9991-41a3-bd69-f4e48df38de0",
   "metadata": {},
   "source": [
    "# Check the observation files now\r\n",
    "Notice that they now have a value for each time step recorded. Not exactly thrilling stuff but lets change our model to vary the pumping and then maybe this will be more interesting to look at.\r\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e94a9ef6-e9cf-436f-a328-8f3a84f0b107",
   "metadata": {},
   "source": [
    "# What if I want to pump more in summer when there is less surface water and less in winter?\r\n",
    "There are different ways to approach this we will do it with a time series because we can vary our stresses at the time step level so no need for lots of stress periods. Previous versions of MODFLOW required that you needed a stress period in order to change any values associated with any transient boundary conditions. The time series needs to extend beyond the activity of the boundary activation that depends on it. For example, if I have a boundary condition active for one year then ideally my time series should extend a day before and a day after to be safe. We will go through the options different options associated with time series files in the workshop or you can read more about how they are configured and what the different variables mean in the MF6io.pdf document. We need to rebuild the well package to change the assigned rates to a time series variable name.\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee5e833e-800e-4c40-9689-6f8a535e706e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove the existing wel package from our simulation\n",
    "gwf.remove_package(\"WEL\")\n",
    "# Now we build a new one starting with the stress period data DICTIONARY\n",
    "# We still use the list of lists but now it is stored inside a dictionary \n",
    "# with the key representing the stress period number\n",
    "\n",
    "wel_period_dict = {} # start an empty dictionary\n",
    "wel1 = [(0,9,9), \"my_rates\",\"Mywell\"] # changed the assigned rate to a string that will be used by the time series package adding in a boundary name\n",
    "wel_period_list = [wel1] # remember that the stress_period data is a list of lists\n",
    "wel_period_dict[1] = wel_period_list # note the dictioanry key is 1, this is stress period two (zero base)\n",
    "\n",
    "wel = flopy.mf6.ModflowGwfwel(gwf, \n",
    "                              boundnames=True, \n",
    "                              stress_period_data=wel_period_dict) # Note we changed this to the dictionary \n",
    "# added the bound and built the object but didn't setup an observation package\n",
    "# Lets do it now using an initialization method\n",
    "\n",
    "# setup the observation dictionary the same as before\n",
    "wel_obs = {\n",
    "    \"wel_flux.csv\": [\n",
    "        (\"Well rate\", \"WEL\", \"Mywell\")]\n",
    "}\n",
    "wel.obs.initialize(continuous=wel_obs)\n",
    "# check the model object to make sure we haven't doubled up on boundaries\n",
    "gwf.get_package_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7349c78d-2155-4c70-91f9-74c9bbf72f7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# here I am going to use a sinusoidal function to generate a variable time series\n",
    "\n",
    "times = np.linspace(1,366,12) # create my sampling periods approximatley monthly over a 366 day period\n",
    "\n",
    "rates = 10.0 * np.cos(2*np.pi*1/366*times)+10.0 \n",
    "# 10.0 represents my amplitude, 2*pi means I want a full wave, \n",
    "# 1/366 is my frequency so the wave only repeats once every 366 time units\n",
    "# times is passed in to say \"give me the values at each of these sampling times\"\n",
    "# the final 10 is an offset so I don't get the positive to negative swing\n",
    "plt.plot(times,rates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b9ed658-fd69-48fc-9e16-ccd222306c5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we need to provide the list of times in model runtime, this means that times are days\n",
    "# \n",
    "ts_data = list(zip(times,rates)) \n",
    "# initialize first time series\n",
    "wel.ts.initialize(\n",
    "    filename=\"flow.wel.ts\", # here I am telling the constructor what filename to use\n",
    "    timeseries=ts_data, # This will become the time and value columns\n",
    "    sfacrecord=-1.0, # this is a scaling factor. Why am I using negative 1.0 here?\n",
    "    time_series_namerecord='my_rates', # this is the name record that the well package has assigned instead of a rate\n",
    "    interpolation_methodrecord=\"linear\", # there are different interpolation approaches. Read about them in the MF6io.pdf doc\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbd8cb03-c25b-4f30-99f0-ed6bb05b0a41",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write all the files\n",
    "wel.write()\n",
    "wel.obs.write()\n",
    "wel.ts.write()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94590565-e28b-48cf-b0b6-b45504393083",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets look at the well file\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,'flow.wel'))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "018c399a-4724-42a1-9120-001a56d5a6af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets look at the well time-series file\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,'flow.wel.ts'))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0673398-3f18-491b-8cf4-68d21e6c1629",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the model object package list\n",
    "gwf.get_package_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15e92762-7f70-4668-a64f-1b2b2796fbd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets run our model again\n",
    "ws3_mod1_trans(sim,gwf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3152c36d-98f3-4fb9-8550-a64cb1b141e2",
   "metadata": {},
   "source": [
    "Hmmm, looks like the model failed because our time series does not extend beyond the end of the simulation. We'll have to try again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d260409-1be0-4d56-a831-368bb900ba4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove the existing wel package from our simulation\n",
    "gwf.remove_package(\"WEL\")\n",
    "# Now we build a new one starting with the stress period data DICTIONARY\n",
    "# We still use the list of lists but now it is stored inside a dictionary \n",
    "# with the key representing the stress period number\n",
    "\n",
    "wel_period_dict = {} # start an empty dictionary\n",
    "wel1 = [(0,9,9), \"my_rates\",\"Mywell\"] # changed the assigned rate to a string that will be used by the time series package adding in a boundary name\n",
    "wel_period_list = [wel1] # remember that the stress_period data is a list of lists\n",
    "wel_period_dict[1] = wel_period_list # note the dictioanry key is 1, this is stress period two (zero base)\n",
    "\n",
    "wel = flopy.mf6.ModflowGwfwel(gwf, \n",
    "                              boundnames=True, \n",
    "                              stress_period_data=wel_period_dict) # Note we changed this to the dictionary \n",
    "# added the bound and built the object but didn't setup an observation package\n",
    "# Lets do it now using an initialization method\n",
    "\n",
    "# setup the observation dictionary the same as before\n",
    "wel_obs = {\n",
    "    \"wel_flux.csv\": [\n",
    "        (\"Well rate\", \"WEL\", \"Mywell\")]\n",
    "}\n",
    "wel.obs.initialize(continuous=wel_obs)\n",
    "\n",
    "times = np.linspace(1,730,24) # create my sampling periods approximatley monthly over a two years.\n",
    "rates = 10.0 * np.cos(2*np.pi*1/365*times)+10.0\n",
    "# we need to provide the list of times in model runtime, this means that times are days\n",
    "# \n",
    "ts_data = list(zip(times,rates)) \n",
    "# initialize first time series\n",
    "wel.ts.initialize(\n",
    "    filename=\"flow.wel.ts\", # here I am telling the constructor what filename to use\n",
    "    timeseries=ts_data, # This will become the time and value columns\n",
    "    sfacrecord=-1.0, # this is a scaling factor. Why am I using negative 1.0 here?\n",
    "    time_series_namerecord='my_rates', # this is the name record that the well package has assigned instead of a rate\n",
    "    interpolation_methodrecord=\"linear\", # there are different interpolation approaches. Read about them in the MF6io.pdf doc\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cac3516-7f0e-4e41-ada3-0f7b52caa778",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the package list\n",
    "gwf.get_package_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b834763-76f4-407a-bad2-150cfd362892",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets run our model again\n",
    "ws3_mod1_trans(sim,gwf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1618be87-7bfe-4e40-924f-e357ccbac110",
   "metadata": {},
   "source": [
    "# Lets plot up some information from our observation files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a94872f-3400-4301-96f0-bbbb8cd720eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting up your observation flux files\n",
    "# We will use pandas to load the observation csv files\n",
    "chd_df = pd.read_csv(os.path.join(model_f,'chd_flux.csv'))\n",
    "wel_df = pd.read_csv(os.path.join(model_f,'wel_flux.csv'))\n",
    "\n",
    "# make the time column a datetime index based on model start date\n",
    "start = pd.to_datetime(\"2023-12-31T00:00:00\") # where did this string come from?\n",
    "time_series_index = start + pd.to_timedelta(chd_df['time'],unit='days') # this is going to give me a time series index\n",
    "\n",
    "# Change the index\n",
    "chd_df.index = time_series_index\n",
    "wel_df.index = time_series_index\n",
    "\n",
    "# drop the time column because we don't want to plot this\n",
    "chd_df.drop(columns='time', inplace=True)\n",
    "wel_df.drop(columns='time', inplace=True)\n",
    "\n",
    "# the well has a no-data value for the steady-state stress period (MF6 no-data = 3.0E30)\n",
    "# We need to convert this to 0.0 or it will cause problems with plotting\n",
    "wel_df.replace(to_replace=3.0E+30, value=np.NaN, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2f1c96c-f2e9-450a-9d76-29f6ab2b5f01",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Easy way to plot with a pandas dataframe if you want to look at your figures\n",
    "ax = chd_df.plot() # will plot all columns in the dataframe as time series\n",
    "ax2=ax.twinx() # add a second y axis\n",
    "wel_df.plot(ax=ax2,linestyle='--',color='red') # change some of the plot styles\n",
    "ax.set_ylabel(\"Flux $m^{3}d^{-1}$\") # Have you seen the $$ used to format math text before?\n",
    "ax2.set_ylabel('Extraction Rate $m^{3}d^{-1}$\"')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07a8d8b5-285a-408a-a99f-f258fe83bc7b",
   "metadata": {},
   "source": [
    "# Multipliers\r\n",
    "Most boundary packages feature a few special types of auxiliary variables. One is a multiplier. The other is actually reserved for transport simulations. Here we will demonstrate the use of an auxiliary variable as a multiplier. The nice thing about using the auxiliary variables is that they can be assigned a time series too so instead of a specified value as shown in the example we can also change the entries from numbers to time series names and add those time series to the existing time series file or create completely separate time series files for the different auxiliary variables. The MF6io.pdf document highlights which variables in each package may be assigned a time series.\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7ed4ffe-8ec3-4c2b-b34a-e8740e98da0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove the existing wel package from our simulation\n",
    "gwf.remove_package(\"WEL\")\n",
    "# Now we build a new one starting with the stress period data DICTIONARY\n",
    "# We still use the list of lists but now it is stored inside a dictionary \n",
    "# with the key representing the stress period number\n",
    "\n",
    "wel_period_dict = {} # start an empty dictionary\n",
    "wel1 = [(0,9,9), \"my_rates\",1.5,100.0,\"Mywell\"] # added mymult=1.5 and concentration=100\n",
    "wel_period_list = [wel1] # remember that the stress_period data is a list of lists\n",
    "wel_period_dict[1] = wel_period_list # note the dictioanry key is 1, this is stress period two (zero base)\n",
    "\n",
    "wel = flopy.mf6.ModflowGwfwel(gwf, \n",
    "                              boundnames=True,\n",
    "                              save_flows=True, # telling it to save flows to the budget file\n",
    "                              auxiliary=['mymult','concentration'], # here I'm assigning 2 auxiliary variables and naming them\n",
    "                              auxmultname='mymult', # here I am saying that the first auxiliary variable called mymult is a multiplier\n",
    "                              stress_period_data=wel_period_dict) # Note we changed this to the dictionary \n",
    "# added the bound and built the object but didn't setup an observation package\n",
    "# Lets do it now using an initialization method\n",
    "\n",
    "# setup the observation dictionary the same as before\n",
    "wel_obs = {\n",
    "    \"wel_flux.csv\": [\n",
    "        (\"Well rate\", \"WEL\", \"Mywell\")]\n",
    "}\n",
    "wel.obs.initialize(continuous=wel_obs)\n",
    "\n",
    "times = np.linspace(1,730,24) # create my sampling periods approximatley monthly over a two years.\n",
    "rates = 10.0 * np.cos(2*np.pi*1/365*times)+10.0\n",
    "# we need to provide the list of times in model runtime, this means that times are days\n",
    "# \n",
    "ts_data = list(zip(times,rates)) \n",
    "# initialize first time series\n",
    "wel.ts.initialize(\n",
    "    filename=\"flow.wel.ts\", # here I am telling the constructor what filename to use\n",
    "    timeseries=ts_data, # This will become the time and value columns\n",
    "    sfacrecord=-1.0, # this is a scaling factor. Why am I using negative 1.0 here?\n",
    "    time_series_namerecord='my_rates', # this is the name record that the well package has assigned instead of a rate\n",
    "    interpolation_methodrecord=\"linear\", # there are different interpolation approaches. Read about them in the MF6io.pdf doc\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64fa13e8-d5da-469f-b5eb-7d33cdcd99d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the package list\n",
    "gwf.get_package_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65b789bb-fbf2-44c4-8e41-7528e41b5923",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets run our model again\n",
    "ws3_mod1_trans(sim,gwf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbb60362-4f64-40b9-a590-ee11cb1da3e3",
   "metadata": {},
   "source": [
    "# What did the multiplier do?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcad4836-7bca-4f35-8c76-340d3170e3d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting up your observation flux files\n",
    "# We will use pandas to load the observation csv files\n",
    "chd_df = pd.read_csv(os.path.join(model_f,'chd_flux.csv'))\n",
    "wel_df = pd.read_csv(os.path.join(model_f,'wel_flux.csv'))\n",
    "\n",
    "# make the time column a datetime index based on model start date\n",
    "start = pd.to_datetime(\"2023-12-31T00:00:00\") # where did this string come from?\n",
    "time_series_index = start + pd.to_timedelta(chd_df['time'],unit='days') # this is going to give me a time series index\n",
    "\n",
    "# Change the index\n",
    "chd_df.index = time_series_index\n",
    "wel_df.index = time_series_index\n",
    "\n",
    "# drop the time column because we don't want to plot this\n",
    "chd_df.drop(columns='time', inplace=True)\n",
    "wel_df.drop(columns='time', inplace=True)\n",
    "\n",
    "# the well has a no-data value for the steady-state stress period (MF6 no-data = 3.0E30)\n",
    "# We need to convert this to 0.0 or it will cause problems with plotting\n",
    "wel_df.replace(to_replace=3.0E+30, value=np.NaN, inplace=True)\n",
    "\n",
    "# Easy way to plot with a pandas dataframe if you want to look at your figures\n",
    "ax = chd_df.plot() # will plot all columns in the dataframe as time series\n",
    "ax2=ax.twinx() # add a second y axis\n",
    "wel_df.plot(ax=ax2,linestyle='--',color='red') # change some of the plot styles\n",
    "ax.set_ylabel(\"Flux $m^{3}d^{-1}$\") # Have you seen the $$ used to format math text before?\n",
    "ax2.set_ylabel('Extraction Rate $m^{3}d^{-1}$\"')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "484676a7-a0c5-458d-9255-9b361c1b2b2e",
   "metadata": {},
   "source": [
    "# What do these methods cover?\n",
    "Using the same approach demonstrated in this workshop you should be able to create all of the following boundary conditions:\n",
    "1. River \n",
    "2. General Head Boundary\n",
    "3. Drain\n",
    "4. Horizontal Flow Barrier\n",
    "5. Well\n",
    "6. Constant head\n",
    "7. Recharge (list based)\n",
    "8. Evapotranspiration (list based)\n",
    "\n",
    "The only difference is the information that needs to be supplied in the stress period data dictionary. For example with the GHB package you would be specifying: \n",
    "\n",
    "a model cell, an hydraulic head, a conductance, auxiliary variables, boundary name\n",
    "\n",
    "There isn't much more to it than that. Selecting your cells can be problematic but we can leverage the intersection objects for that."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7291dcd7-0978-4c16-a194-20b431129834",
   "metadata": {},
   "source": [
    "# Selcting Cells with grid intersection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6577c406-0f0d-41b1-ae50-fe47372839f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Selecting cells using Flopy intersection methods\n",
    "from shapely.geometry import (\n",
    "    LineString,\n",
    "    Point,\n",
    "    Polygon,\n",
    ")\n",
    "from flopy.utils import GridIntersect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9abdfae-fe1b-4c34-ab78-2cfb558ced0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# examples of shapefile geometries obtained from shapley or geopandas\n",
    "# you can create these shapefiles in a GIS if you have a projected grid\n",
    "# if there is sufficient time then we will demonstarte this too.\n",
    "poly = Polygon(shell=[(3, 3),(15.6, 7.7),(14.5, 12.9),(4.7, 10.2)])\n",
    "line = LineString(coordinates=[(5,4),(19,18)])\n",
    "point = Point([17.6,5.7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a91cea83-5e7a-4207-9079-598125e27144",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create an intersect object\n",
    "ix = GridIntersect(mg, method=\"structured\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "794488cd-d498-44f7-beb1-6dd83f9b5ca8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#perform all the intersections\n",
    "result1 = ix.intersect(poly)\n",
    "result2 = ix.intersect(line)\n",
    "result3 = ix.intersect(point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc30538b-dd67-4f38-8565-834ffd671ce5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the cell ids covered by the polygon\n",
    "result1.cellids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d5bbbc2-7010-451b-8bbc-9cc753ae719a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use them to make a recharge boundary condition\n",
    "# annual rainfall mean = 1220 mm\n",
    "# we will start with 1% of annual rainfall\n",
    "# estimate of 0.5% to 2% of annual\n",
    "rain_rate = 1220/365/1000 # (m/d)\n",
    "min_rate = 0.005*rain_rate\n",
    "max_rate = 0.03*rain_rate\n",
    "# our preferred value which is central to log of min and max\n",
    "pv =  10**(np.log10(min_rate)+((np.log10(max_rate)-np.log10(min_rate))/2)) \n",
    "\n",
    "rch_list = [[(0,x[0],x[1]),pv,'poly_rch'] for x in result1.cellids] # make the list with the cell ids. \n",
    "#Why is l,r,c assigned like this (0,x[0],x[1])?\n",
    "rch_period = {} # start the dictionary\n",
    "rch_period[0]=rch_list \n",
    "rch = flopy.mf6.ModflowGwfrch(\n",
    "    gwf,\n",
    "    boundnames=True,\n",
    "    save_flows=True,\n",
    "    stress_period_data=rch_period,\n",
    ")\n",
    "# setup the observation dictionary the same as before\n",
    "rch_obs = {\n",
    "    \"rch_flux.csv\": [\n",
    "        (\"Recharge\", \"RCH\", \"poly_rch\")]\n",
    "}\n",
    "rch.obs.initialize(continuous=rch_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf5c9357-f35d-4969-a7c7-070eaed02fd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets get a modelgrid object and have a look at the grid\n",
    "mg = gwf.modelgrid\n",
    "\n",
    "# create a figure instance to save to our plots folder\n",
    "fig,ax = plt.subplots(figsize=(7,5)) # we are creating a figure object here so that we can dictate size note there are mutiple ways to do this\n",
    "\n",
    "# Create the mapview object of the modelgrid and tell it to use the axes we created with our figure\n",
    "pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) # note by not specifying a layer here it will assume layer 1\n",
    "\n",
    "# say that you want the Mapview object to plot the grid with linewidth = 0.3 and colour the lines black\n",
    "pmv.plot_grid(ax=ax, lw=0.3, color=\"black\")\n",
    "\n",
    "# Now I also want to plot the boundary condition so I specify the package object to plot\n",
    "pmv.plot_bc(package=rch, color='blue')\n",
    "\n",
    "\n",
    "ax.set_title('RCH Boundary Cells')\n",
    "ax.set_xlabel('Eastings')\n",
    "ax.set_ylabel('Northings')\n",
    "# lets save it to our plots folder\n",
    "figname = os.path.join(plots_f,'RCH_bounadry_cells_LR.png') # note use of the plots folder path.\n",
    "# If you want to change the file format then change the extension from .png to pdf or just do both\n",
    "fig.savefig(figname,dpi=300)\n",
    "figname = os.path.join(plots_f,'RCH_bounadry_cells_LR.pdf')\n",
    "fig.savefig(figname,dpi=300) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36354d01-4e9e-4930-ba9d-3f538d1450b1",
   "metadata": {},
   "source": [
    "# Depending on time.\n",
    "If there is time available then we will have a look at a projected grid with inersection and setting up complex transient boundary conditions that would take a long time with GUI's but is relativley quick with a script. If no time is available then we will leave it for later in the series."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d4c74ad-4923-4457-b4e7-2fa7dc1c5142",
   "metadata": {},
   "source": [
    "# Thats all folks!\n",
    "Next session will be on configuring array style properties and boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "984ad090-60a3-4e16-b711-9c1a0750d3b1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}