Workshop_4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f6b92ba0-7549-4a25-9b60-06d182512f2e",
   "metadata": {},
   "source": [
    "# Workshop 4\n",
    "Hi and welcome to the fourth workshop. In this session we will be focusing on creating input packages that use array style input. We've covered some array style input previously when building our grids so hopefully this will be a nice extension. Two boundary conditions, the ET and recharge package allow for both list and array cofiguration. Because MF6 allows for multiple instances of the same type of stress package to be included in a model you can have both an array style recharge package combined with a list one if you need that functionality. The different property packages also use array style input. These include hydraulic conductivity with the NPF package, storage with the STO package, skeletal storage, compaction and subsidence with the CSUB package, the initial conditions IC package and many groundwater transport packages use array style input. Many array style input packages allow you to forego providing an array if values are consistent across the entire array. Thankfully, once you've setup a few of these packages the appraoch will become familiar and you should be able to use the same methods for all. At the end of the session we will cover the output control package and how to configure it plus the model level observation package (as opposed to the stress level observations we covered in the last session)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe7c3abd-b61c-4853-a79f-4ac1b4fe6f8c",
   "metadata": {},
   "source": [
    "# Imports\n",
    "These should need no explanation by now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6a96aab-c3cd-4e74-b86e-fc47bcdfad44",
   "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",
    "from flopy.discretization import StructuredGrid, VertexGrid\n",
    "\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": "216f3e3c-5eba-4e82-b4b7-d4128a5160ab",
   "metadata": {},
   "source": [
    "# Folder setup\n",
    "Again nothing here you haven't seen before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "075afe21-7100-4a6c-96c1-f1ba2b5dbc35",
   "metadata": {},
   "outputs": [],
   "source": [
    "ws4 = os.path.join('workshop_4') # here we are making a path not creating the folder\n",
    "gis_f = os.path.join(ws4,'GIS') # creating a sub-directory path for our gis input/output\n",
    "model_f = os.path.join(ws4,'model') # creating a sub-directory path for our model input/output\n",
    "plots_f = os.path.join(ws4,'plots') # creating a sub-directory path for our plots\n",
    "for path in [ws4,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": "3dbb1ab5-c42b-4445-aa38-469a814f7387",
   "metadata": {},
   "source": [
    "# Build a transient model\r\n",
    "Last session we started off with a steady state model. This time we will jump straight to a transient model. Keep an eye out for a few other options being activated in the simulation level objects that you may not have seen before\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7402b0d-9030-4e95-9acc-3301fb370d41",
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup sim\n",
    "sim_name = \"MySim\" \n",
    "sim = flopy.mf6.MFSimulation(sim_name=sim_name, \n",
    "                             exe_name=\"mf6\",\n",
    "                             verbosity_level=1,\n",
    "                             sim_ws=model_f) \n",
    "\n",
    "# build timing data (SS plus 10 years monthly transient stress period)\n",
    "start_date = \"2023-12-31\" # this will be our SS date\n",
    "dates = pd.date_range('2024-01-01','2034-01-01', freq='MS').tolist() # note this will be our transient period dates\n",
    "perlens = [(dates[x]-dates[x-1]).days for x in range(1,len(dates))]\n",
    "stp = 1 \n",
    "_ = [(x,stp,1) for x in perlens]\n",
    "pdata = [(1,1,1), *_, (36525, 1200, 1)]\n",
    "perlens = [x[0] for x in pdata]\n",
    "numper = len(pdata)\n",
    "\n",
    "# build a stress period timing dataframe\n",
    "dates = [start_date,*dates] \n",
    "df = pd.DataFrame() \n",
    "df['Date'] = dates \n",
    "df['SP'] = range(1,len(dates)+1) \n",
    "df['Flopy_SP'] = range(len(dates)) \n",
    "df['Incremental'] = perlens\n",
    "df['Cumulative'] = np.cumsum(perlens) \n",
    "df.to_csv(os.path.join(model_f,'model_timing.csv'),index=None) \n",
    "modtime_df = df.copy()\n",
    "\n",
    "# setup tdis\n",
    "tdis = flopy.mf6.ModflowTdis(sim,\n",
    "                             time_units='days',\n",
    "                             nper=numper,\n",
    "                             perioddata=pdata,\n",
    "                             start_date_time=start_date) \n",
    "\n",
    "# IMS\n",
    "ims = flopy.mf6.ModflowIms(sim, complexity='MODERATE', \n",
    "                           csv_inner_output_filerecord='inner.csv', \n",
    "                           csv_outer_output_filerecord='outer.csv', \n",
    "                           outer_maximum=500, \n",
    "                           inner_maximum=500, \n",
    "                           outer_dvclose=0.01, \n",
    "                           inner_dvclose=0.001) \n",
    "\n",
    "# GWF\n",
    "model_name = 'flow' \n",
    "gwf = flopy.mf6.ModflowGwf(sim, \n",
    "                           modelname=model_name, \n",
    "                           save_flows=True, \n",
    "                           newtonoptions=\"under_relaxation\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32b3de38-0313-4e53-b792-78945beaf838",
   "metadata": {},
   "outputs": [],
   "source": [
    "# copy the shapefiles across to our project working directory for GIS\n",
    "# again this should be familiar\n",
    "shp_path = os.path.join('files','disv_shapefiles') # path to shapefiles for this example\n",
    "flist = [x for x in os.listdir(shp_path)] # create a list of all the shapefiels\n",
    "for file in flist:\n",
    "    shutil.copyfile(os.path.join(shp_path,file),os.path.join(gis_f,file)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1060d8b-10b6-4e30-a842-6d29941e39aa",
   "metadata": {},
   "source": [
    "# Triangle\r\n",
    "Below we make a similar grid to what we used previously with Gridgen only this time we use Triangle and build a Voronoi grid instead of a quadtree refined grid I've commented it out (select all then hit ctrl+/ to block-comment or reverse) it out for now but felt it was worthwhile including it. \r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "332c3fb5-c0f4-4c6b-87f6-2bb011368493",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # First we need to import some libraries and objects\n",
    "# from flopy.utils.triangle import Triangle as Triangle\n",
    "# from flopy.utils.voronoi import VoronoiGrid\n",
    "# import geopandas as gpd\n",
    "\n",
    "# # get information for model domain * must be the first polygon added to triangle object \n",
    "# # the first polygon that is added to Triangle will be the domain this is a requirement of the mesh creation\n",
    "# ad = os.path.join(gis_f,\"model_bounds_poly.shp\") # get the path to the model boundary shape\n",
    "# ad_df = gpd.read_file(ad) # use geopandas to read the shapefile into a geo-dataframe *** this is new\n",
    "# # We need to pass a shapely polygon object to Triangle and also a point inside the polygon for limiting the cell area\n",
    "# ad_shape = ad_df.geometry[0] # access the shapley geometry object from the geopandas dataframe\n",
    "# ad_point = (ad_df.centroid.x[0],ad_df.centroid.y[0]) #get a point in the center of the polgon as a tuple of x\n",
    "# ad_area = 1280**2 # setting the maximum cell area in the model to be comparable with Gridgen Note this is the maximum allowed meaning anything less goes.\n",
    "\n",
    "# # You can also specify node positions directly. This is conveneint for boundary condition assignment but requires some GIS work upfront\n",
    "# # So here we get explicit nodes that we created in GIS.\n",
    "# mc = os.path.join(gis_f,\"channel_points.shp\") # this is a point shapefile that I use to specify node coordinates along my streams \n",
    "# mc_df = gpd.read_file(mc) # use geopandas to read the point shapefile\n",
    "# mc_point_array = np.array(list(zip(mc_df.geometry.x,mc_df.geometry.y))) # create an array of x,y tuples to provide to Triangle, think of this a column of (x,y)\n",
    "\n",
    "\n",
    "# # This is the same process used for the model polygon but will now be for a refined area around the project\n",
    "# vb = os.path.join(gis_f,\"vgrid_buffer.shp\")\n",
    "# vb_df = gpd.read_file(vb)\n",
    "# vb_shape = vb_df.geometry[0]\n",
    "# vb_point = (vb_df.centroid.x[0],vb_df.centroid.y[0])\n",
    "# vb_area = (1280**2)/(4*7) # note here I'm limiting the cell area to be comparable to 6 levels of equivalent quadtree refinement\n",
    "\n",
    "\n",
    "# # Now we start the Triangle object, spot where all the information we extracted above is used\n",
    "# tri = Triangle(angle=30,nodes=mc_point_array,model_ws=model_f) # What does angle=30 mean?\n",
    "\n",
    "# # then add in our polygons, note which one apears first\n",
    "# tri.add_polygon(ad_shape)\n",
    "# tri.add_polygon(vb_shape)\n",
    "# # then we mark each polygon with a point inside them which declares them as a region\n",
    "# # and we limit the cell surface area in that region\n",
    "# tri.add_region(ad_point, 0, maximum_area=ad_area)\n",
    "# tri.add_region(vb_point, 1, maximum_area=vb_area)\n",
    "# tri.build() # the build process overall is not that disimilar to gridgen\n",
    "\n",
    "# vor = VoronoiGrid(tri) # create a voronoigrid object from the Triangle one\n",
    "# gridprops = vor.get_gridprops_vertexgrid() # recall the difference between modelgrid properties and disv properties\n",
    "# idomain_vor = np.ones((1, vor.ncpl), dtype=int) # this is a dummy idomain for now, we will assigne a three layer one later\n",
    "# voronoi_grid = VertexGrid(**gridprops, nlay=1, idomain=idomain_vor) # this is akin to the modelgrid is not yet registered to a model object\n",
    "\n",
    "# # make a figure\n",
    "# fig = plt.figure(figsize=(20, 20))\n",
    "# ax = plt.subplot(1, 1, 1, aspect=\"equal\")\n",
    "# voronoi_grid.plot(ax=ax, facecolor=\"none\")\n",
    "# ax.ticklabel_format(style='plain') # test what happens if you don't use this switch\n",
    "# ax.set_title('DISV Voronoi 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,'DISV_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,'DISV_model_grid.pdf')\n",
    "# fig.savefig(figname,dpi=300) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74d15dc7-bdf3-4ac8-8f4c-8dc74de14197",
   "metadata": {},
   "source": [
    "# Our Gridgen grid from Workshop 2\n",
    "The following cells were copied from the previous workshop so need no real explanation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36528e3a-cb0c-4c28-a94c-02203c2e3fca",
   "metadata": {},
   "outputs": [],
   "source": [
    "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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "979a33d4-555e-4dbd-9ea3-f61638113fa4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# start with our basic structured grid\n",
    "nlay = 3\n",
    "nrow = 34\n",
    "ncol = 44\n",
    "delr = delc = 1280.0\n",
    "botm = np.zeros((nlay, nrow, ncol), dtype=np.float32)\n",
    "top = np.zeros((1, nrow, ncol), dtype=np.float32)\n",
    "idom = np.ones((nlay, nrow, ncol), dtype=np.float32)\n",
    "botm[0, :, :] = 390.0\n",
    "botm[1,:,:] = 380.0\n",
    "botm[2,:,:] = -170.0\n",
    "top[0,:,:] = 460.0\n",
    "\n",
    "\n",
    "# Note we start with a structured DIS grid despite aiming for a DISV grid.\n",
    "dis = flopy.mf6.ModflowGwfdis(\n",
    "    gwf,\n",
    "    nlay=nlay,\n",
    "    nrow=nrow,\n",
    "    ncol=ncol,\n",
    "    delr=delr,\n",
    "    delc=delc,\n",
    "    top=top,\n",
    "    botm=botm,\n",
    "    xorigin=729425,\n",
    "    yorigin=947000,\n",
    "    length_units='meters',\n",
    "    angrot=0,\n",
    "    idomain = idom\n",
    ")\n",
    "dis.export(os.path.join(gis_f,'disv.shp'))\n",
    "# check that the exported disv shapefile covers the polygon for the model boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90903f3c-0b5b-4081-9fba-5efba997e1fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import Polygon\n",
    "g = Gridgen(dis)\n",
    "dam = os.path.join(gis_f,\"dam_buffer\")\n",
    "chanel = os.path.join(gis_f,\"my_channels\")\n",
    "tsf = os.path.join(gis_f,\"tsf_buffer\")\n",
    "wels = os.path.join(gis_f,\"Wells_buffered\")\n",
    "pit1500 = os.path.join(gis_f,\"pits_buffer_1500\")\n",
    "pit1000 = os.path.join(gis_f,\"pits_buffer_1000\")\n",
    "pit500 = os.path.join(gis_f,\"pits_buffer_500\")\n",
    "mod_bnd = os.path.join(gis_f,\"model_bounds\")\n",
    "act_dom = os.path.join(gis_f,\"model_bounds_poly\")\n",
    "\n",
    "g.add_refinement_features(chanel, \"line\", 3, layers=[0,1,2])\n",
    "g.add_refinement_features(wels, \"polygon\", 3, layers=[0,1,2])\n",
    "g.add_refinement_features(dam, \"polygon\", 3, layers=[0,1,2])\n",
    "g.add_refinement_features(tsf, \"polygon\", 3, layers=[0,1,2])\n",
    "g.add_refinement_features(pit1500, \"polygon\", 3, layers=[0,1,2])\n",
    "g.add_refinement_features(mod_bnd, \"line\", 3, layers=[0,1,2])\n",
    "g.add_refinement_features(pit1000, \"polygon\", 4, layers=[0,1,2])\n",
    "g.add_refinement_features(pit500, \"polygon\", 5, layers=[0,1,2])\n",
    "g.add_active_domain(act_dom,layers=[0,1,2])\n",
    "g.build()\n",
    "\n",
    "grd_files = [file for file in os.listdir('.') if file.startswith(\"qtgrid\")]\n",
    "for file in grd_files:\n",
    "    shutil.copyfile(file,os.path.join(gis_f,file))\n",
    "\n",
    "gridprops_vg = g.get_gridprops_vertexgrid()\n",
    "vgrid = flopy.discretization.VertexGrid(**gridprops_vg)\n",
    "fig,ax = plt.subplots(figsize=(12,12))\n",
    "vgrid.plot(ax=ax)\n",
    "ax.set_ylabel('Northing')\n",
    "plt.title('Model Grid')\n",
    "figname = os.path.join(plots_f,'model_grid.png') \n",
    "fig.savefig(figname,dpi=300)\n",
    "figname = os.path.join(plots_f,'model_grid.pdf')\n",
    "fig.savefig(figname,dpi=300) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6438549-d941-4928-9988-9ae5a38a19c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now lets get the properties for a disv object. Note the method is different and is specific to a disv type object\n",
    "gridprops_disv = g.get_gridprops_disv()\n",
    "\n",
    "# rebuild gwf *** IF YOU DON'T REBUILD GWF IT WILL FAIL *** \n",
    "gwf = flopy.mf6.ModflowGwf(sim, modelname=model_name, save_flows=True, newtonoptions=\"under_relaxation\")\n",
    "\n",
    "# Note use of mf6.ModflowGwfdisv, passing in the model object (gwf) and unpacking the grid properties dictionary\n",
    "disv = flopy.mf6.ModflowGwfdisv(gwf,angrot=0,length_units=\"METERS\", **gridprops_vg)\n",
    "# ignore the warning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "932f6048-889d-4720-b1c1-cebe82c33136",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a path to where we want our new folder\n",
    "grid_path = os.path.join(ws4,'Gridgen')\n",
    "# now make the folder by checking if it exists first, deleting it, then making it\n",
    "if os.path.exists(grid_path): # here we are asking if the path exists on the computer. \n",
    "    shutil.rmtree(grid_path)# if it does exist, delete it and all the files in it plus any subdirectories\n",
    "    os.mkdir(grid_path) # then remake it\n",
    "else:\n",
    "    os.mkdir(grid_path) # if it doesn't exist then make the folder\n",
    "# now make a list of all the files we want to move\n",
    "# We will use the start of the different filenames to identify them\n",
    "flist = [] # start with an empty list\n",
    "for pref in ['qtg', 'quadtree', '_gridgen',]: # kick off a for loop with a list of file prefixes\n",
    "    temp_list = [x for x in os.listdir() if x.startswith(pref)] # make a temporary list of the files that start with this loops prefix\n",
    "    flist = [*flist,*temp_list] # unpack the existing list of files and the temporary list of files into the file list\n",
    "# Before we move the files you should check that you haven't included any files you don't want to move.\n",
    "flist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c782846-f013-4790-80df-7531466793bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# If that looks good then we move the files to our new folder\n",
    "for file in flist:\n",
    "    shutil.move(file,grid_path) # note use of move as opposed to copyfile which we used earlier. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48364473-540b-4a7a-a4f1-647971600819",
   "metadata": {},
   "source": [
    "# Methods to build arrays for your model\n",
    "Since this workshop will focus on building array style boundary conditions and property packages then it may be useful to examine a few methods for building arrays for your model to begin with. We will start with some basics which we already demonstrated in previous workshops. The modelgrid object gives you access to the information you need to size your arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6aa6bc73-56fd-4177-8404-6a293dbcf8ed",
   "metadata": {},
   "source": [
    "# Basic arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d659239-806b-4702-b4d9-99a5dc488382",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting up basic arrays\n",
    "# make sure your modelgrid is updated\n",
    "mg = gwf.modelgrid\n",
    "\n",
    "# Start with the size of the array you need\n",
    "l1_array = np.ones_like(mg.top) # makes an array filled with 1.0 having the size for layer 1 often used for recharge and ET array creation\n",
    "# l1_array = np.zeros_like(mg.top) 3\n",
    "\n",
    "mod_array = np.ones_like(mg.botm) # makes an array with the size for the entire model often used for property assignments\n",
    "\n",
    "# check shapes of arrays\n",
    "print(np.shape(l1_array), np.shape(mod_array))\n",
    "\n",
    "# But what if I want to assign different values to a specific layer\n",
    "# Easy, start with a list of the values you want to have in each layer\n",
    "layer_props = [0.5,0.05,0.005] \n",
    "# then initialize an array with ones \n",
    "my_array = np.ones_like(mg.botm)\n",
    "# then loop through each layer\n",
    "for i,j in enumerate(layer_props): # the first array \n",
    "    my_array[i]=j\n",
    "my_array\n",
    "\n",
    "# Thats great and all, but what if I also want to have different values within the same layer."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11a3ae91-9711-468d-bb1a-1b3b08c556bf",
   "metadata": {},
   "source": [
    "# A zone within a layer\r\n",
    "You already have the array but you want to assign a values to specific cells in the layer. You can do this with shapefiles via an intersection object with the modelgrid like we demonstrated in the last workshop. The cell selection can be used to index the array. All you have to be aware of is how this changes depending on the grid type you are dealing with. This was demonstrated previously with a DIS grid and here we will use a DISV grid\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e464a796-d6fa-40c5-b651-a5d3dd56cba9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from flopy.utils import GridIntersect\n",
    "import geopandas as gpd\n",
    "# create an intersect object\n",
    "ix = GridIntersect(mg, method=\"vertex\")\n",
    "my_shape = os.path.join(gis_f,\"pits_buffer_1500.shp\")\n",
    "my_poly = gpd.read_file(my_shape).geometry[0]\n",
    "my_cells = ix.intersect(my_poly)\n",
    "my_node_idx = [x[0] for x in my_cells] # to get just the node indices as a list\n",
    "\n",
    "# the index of the cells matches the index in the array so you can edit it\n",
    "# note this works with DIS and DISV but care needs to be taken when working with DISU so always check by plotting\n",
    "my_array[0][my_node_idx] = 1.0 # changing the nodes in layer 1 inside the \"zone\" to have a value of 1.0\n",
    "\n",
    "# plot to check\n",
    "mg = gwf.modelgrid\n",
    "fig,ax = plt.subplots(figsize=(8,8)) # we are creating a figure object here so that we can dictate size note there are mutiple ways to do this\n",
    "pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) # note by not specifying a layer here it will assume layer 1\n",
    "pc = pmv.plot_array(my_array[0])\n",
    "ax.set_aspect('equal')\n",
    "ax.set_title('My Zone')\n",
    "ax.set_xlabel('Eastings')\n",
    "ax.set_ylabel('Northings')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1eb13aae-dcaa-46a8-afcd-dfcdb3e753b0",
   "metadata": {},
   "source": [
    "# Raster sampling to the modelgrid\r\n",
    "We did this previously in Workshop 2 but it demonstrates a method that can be used to build a layer array for any property of boundary condition. It should be clear how this method can be adopted to assigning values on a layer-by-layer basis using rasters of various properties. In the example below a DEM raster is used but this could easily be substituted with rooting depth approximation derived from vegetation mapping to inform your ET package\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afa26799-be64-442a-b9f1-1a2782860a0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we need to use the grid information to make our layers\n",
    "# first we start with the top of the model\n",
    "# this is identical to what was done in Workshop 2 only now with a voronoi grid\n",
    "from flopy.utils import Raster\n",
    "\n",
    "topo_fyl = os.path.join('.','files','filled_dem.tif')\n",
    "\n",
    "rio1 = Raster.load(topo_fyl)\n",
    "fig = plt.figure(figsize=(8, 5))\n",
    "ax = fig.add_subplot(1, 1, 1, aspect=\"equal\")\n",
    "\n",
    "mg=gwf.modelgrid\n",
    "\n",
    "ax = rio1.plot(ax=ax)\n",
    "plt.colorbar(ax.images[0], aspect=30)\n",
    "pmv = flopy.plot.PlotMapView(modelgrid=mg)\n",
    "pmv.plot_grid(ax=ax, lw=0.3, color=\"black\")\n",
    "top_data = rio1.resample_to_grid(\n",
    "    mg, band=rio1.bands[0], method=\"nearest\"\n",
    ")\n",
    "\n",
    "# slice elevation peaks because they won't apply to GW anyway\n",
    "top_data[top_data>450.0]=450.0\n",
    "\n",
    "# Check the size. Note it is for 1 layer only despite the modelgrid having mutiple layers\n",
    "print(np.shape(top_data))\n",
    "\n",
    "ax.set_title('Topography (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\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9846db3-7f94-4d25-a77b-929e7fc69b45",
   "metadata": {},
   "source": [
    "# Scale mapping to build arrays\r\n",
    "The approach demonstrated below can be useful in situations where you have sufficient justification to assume a scaled correlation to exist between different properties. In the example below we build an array of bottom elevations using this approach. Essentially, what we are saying here is that we know what the approximate maximum and minimum thicknesses are of the different aquifers and where the er is greater elevation we expect increased thickness in the top two aquifers. We can also use this approach to assume a scaled correlation to exist between hydraulic conductivity and storage\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0aa7cc9-433d-4617-9695-22addeca7bb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using a range mapping function to fudge some thicknesses \n",
    "\n",
    "# creating a function to apply to each model cell it \n",
    "# basically you provide the max and min for a known range (mx1, mn1)\n",
    "# and the new range (mx2, mn2) you want to map a value from the known range (x)\n",
    "# effectivley you are assuming a direct correlation\n",
    "def scale_me(mx1,mn1,mx2,mn2,x):\n",
    "    r1 = mx1-mn1\n",
    "    r2 = mx2-mn2\n",
    "    return((((x-mn1)*r2)/r1)+mn2)\n",
    "vf = np.vectorize(scale_me) # We vectorize the function just to make it quicker\n",
    "\n",
    "# now I'm going to map the range of known surface elevations to an\n",
    "# approximated range of thickness that I have some idea exists in this region \n",
    "# effectivley what I am saying is that where there is greater elevationI have greater thickness\n",
    "tmax = np.max(top_data)\n",
    "tmin = np.min(top_data)\n",
    "l1max = 60.0\n",
    "l1min = 30.0\n",
    "l1range = l1max-l1min\n",
    "l1_thickness = vf(tmax,tmin,l1max,l1min,top_data) # this is an array of thickness for layer 1 directly correlated with elevation\n",
    "\n",
    "l2max = 65.0\n",
    "l2min = 50.0\n",
    "l2_thickness = vf(tmax,tmin,l2max,l2min,top_data) # this is an array of thickness for layer 2 also directly correlated with elevation\n",
    "\n",
    "# need to build bottom elevations for the whole model as an array\n",
    "new_botms = np.ones_like(mg.botm) # create an array for the bottoms\n",
    "new_botms[0] = top_data - l1_thickness # recall that top_data has array size of 1 layer\n",
    "new_botms[1] = new_botms[0] - l2_thickness\n",
    "new_botms[2] = new_botms[1]-370.0 \n",
    "l3_thickness = new_botms[1] - new_botms[2] # this is an array of thickness for layer 3 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "813c201d-71e0-44f0-9964-5766eaf8a460",
   "metadata": {},
   "outputs": [],
   "source": [
    "# update disv so that it includes the new elevations\n",
    "disv.botm = new_botms\n",
    "disv.top = top_data\n",
    "#update modelgrid the model grid whenever we update dis objects\n",
    "mg = gwf.modelgrid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beef2d04-88ea-496c-9dbd-65fadb59ab30",
   "metadata": {},
   "source": [
    "# Plotting layer information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "396806f9-e00e-40b3-b638-e44c0f1f186c",
   "metadata": {},
   "outputs": [],
   "source": [
    "min = np.min(new_botms)\n",
    "max = np.max(new_botms)\n",
    "\n",
    "for i in range(nlay):\n",
    "    fig = plt.figure(figsize=(8, 5))\n",
    "    ax = fig.add_subplot(1, 1, 1, aspect=\"equal\")\n",
    "    pmv = flopy.plot.PlotMapView(modelgrid=mg)\n",
    "    pc = pmv.plot_array(new_botms[i], masked_values=[1e+30], vmin=min, vmax=max)\n",
    "    cbar = plt.colorbar(pc, aspect = 30)\n",
    "    ax.set_title(f'Base of Layer {i+1} (mAMSL)')\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()\n",
    "    fpath = os.path.join(plots_f,f\"Bottom_layer{i+1}.png\")\n",
    "    fig.savefig(fpath,dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad190df1-3d38-40d7-b629-0d492be333a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "min = np.min(mg.cell_thickness)\n",
    "max = np.max(mg.cell_thickness)\n",
    "for i in range(nlay):\n",
    "    fig = plt.figure(figsize=(8, 5))\n",
    "    ax = fig.add_subplot(1, 1, 1, aspect=\"equal\")\n",
    "    pmv = flopy.plot.PlotMapView(modelgrid=mg)\n",
    "    pc = pmv.plot_array(mg.cell_thickness[i], masked_values=[1e+30], vmin=min, vmax=max)\n",
    "    cbar = plt.colorbar(pc, aspect = 30)\n",
    "    ax.set_title(f'Thicknes of Layer {i+1} (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\n",
    "    plt.tight_layout()\n",
    "    fpath = os.path.join(plots_f,f\"Thickness_layer{i+1}.png\")\n",
    "    fig.savefig(fpath,dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a833b0dc-bb94-44a6-b36a-6e3868b3e9ba",
   "metadata": {},
   "source": [
    "# The NPF package for hydraulic conductivity\r\n",
    "The NPF package is used to assign hydraulic conductivity values to the model. This includes both horizontal and vertical. There are many options that you can set in the NPF package which may very large effects on your solutions so it is worthwhile reading what all the options do for you. At a bare minimum you need to know what \"k\", \"k33\" and \"icelltype\" are. That said configuring the package through Flopy is relatively straight forward. Options are generally set using either True, None/False or a string. The data is provided as an array or alternatively as constants for whole arrays or layers. Note unlike boundary conditions only one NPF package can be assigned to a model\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47b26a0b-d8cd-4eea-834d-46f436bebc6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign fixed values to whole model\n",
    "\n",
    "npf = flopy.mf6.ModflowGwfnpf(gwf,k=1.0,k33=0.1,icelltype=1)\n",
    "npf.write()\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,f'{model_name}.npf'))]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "514ea793-41c4-4f1b-bf7a-9f6ac3046590",
   "metadata": {},
   "source": [
    "Now we can pass in values for a complete layer as a list with a value for each layer. This will invoke the \"layered\" option for input file construction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91fdbd37-c6f6-4799-876c-b1996c5c71c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now with values for specific layers\n",
    "k_layer_list = [1.0, 0.1, 0.01]\n",
    "k33_layer_list = [x/10 for x in k_layer_list]\n",
    "icell_list = [1,-1,0]\n",
    "npf = flopy.mf6.ModflowGwfnpf(gwf,\n",
    "                              k=k_layer_list,\n",
    "                              k33=k33_layer_list,\n",
    "                              icelltype=icell_list,\n",
    "                              thickstrt=True)\n",
    "npf.write()\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,f'{model_name}.npf'))]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c1ee75c-bdef-4ede-84f6-e1bcf386f1b4",
   "metadata": {},
   "source": [
    "Or alternatively you can build the arrays using any of the methods shown previously and pass them in explicitly. The input file will change quite a bit if you do this. The following achieves the exact same outcome as the previous method but passes in the arrays instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6a99e2e-ed13-4a6d-affc-609dc8f176a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "kx_array = np.ones_like(mg.botm)\n",
    "kv_array = kx_array.copy()\n",
    "icell_array = kx_array.copy()\n",
    "for i,(j,k,l) in enumerate(zip(k_layer_list,k33_layer_list,icell_list)): \n",
    "    kx_array[i]=j\n",
    "    kv_array[i]=k\n",
    "    icell_array[i]=l\n",
    "\n",
    "npf = flopy.mf6.ModflowGwfnpf(gwf,\n",
    "                              k=kx_array,\n",
    "                              k33=kv_array,\n",
    "                              icelltype=icell_array,\n",
    "                              save_flows=True,\n",
    "                              thickstrt=True,\n",
    "                              xt3doptions='RHS')\n",
    "npf.write()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "119e40d9-6cc5-40f4-a97e-a18c6fc0d801",
   "metadata": {},
   "source": [
    "# Transient Hydraulic conductivity\r\n",
    "Aquifer properties can be varied during a simulation using the TVK package. This is demonstrated below using the assumption that the hydraulic conductivity needs to vary over time in the same zone that we created earlier but through layers one and two. The package creation is akin to a boundary condition but because you can only have one NPF package you can also only have one TVK package. Here we will use a time series to control how the hydraulic conductivity has to vary over time. MF6 will then adjust the values for the inter-cell conductance calculations (what the hydraulic conductivity arrays are used for)\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0df8b0bb-9e50-43aa-8588-c99438339cbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# just like with our boundary conditions we need to make a stress period dictionary\n",
    "# we start with the data needed for the first\n",
    "pdata=[]\n",
    "for lay in range(2):\n",
    "    for node in my_node_idx:\n",
    "        pdata.append(((lay,node),'k','trans_k')) # 'k' is the property I want to change, 'trans_k' will be the time series that I want to use\n",
    "        pdata.append(((lay,node),'k33','trans_k33')) # 'k33' is the property I want to change, 'trans_k33' will be the time series that I want to use\n",
    "\n",
    "tvk_period = 60 # only want the hydraulic conductivity to change from stress period 60 onwards\n",
    "\n",
    "tvk_pdata = {} # our stress period dictionary\n",
    "for key in range(tvk_period,tvk_period+13): # This means that the package will be active for 1 year.\n",
    "    tvk_pdata[key] = pdata\n",
    "tvk = flopy.mf6.ModflowUtltvk(npf,  #note how we pass in the package here and not the model object or the simulation object\n",
    "                              perioddata=tvk_pdata,\n",
    "                              filename=\"{}.tvk\".format(model_name))\n",
    "\n",
    "# make the time series data. Recall that timeseries require modelruntime values for timing\n",
    "# here we will use the timing dataframe we created earlier specifcally the cumulative column\n",
    "#model_time = np.cumsum(tdis.perioddata.array['perlen'])\n",
    "stime = modtime_df.loc[modtime_df['Flopy_SP']==60,'Cumulative'].iloc[0] # model run time in days at stress period 60\n",
    "sptime = modtime_df['Cumulative'].tolist()[-1] # model run time in days at end of model run\n",
    "ts_data = [(1.0,10.0,10.0),(stime,10.0,10.0),(sptime+1.0,1.0,1.0)] \n",
    "# note you must have this start at t=1.0 and end at the end of the model run\n",
    "\n",
    "\n",
    "# initialize first time series\n",
    "tvk.ts.initialize(\n",
    "    filename=\"tvk.ts\",\n",
    "    timeseries=ts_data,\n",
    "    time_series_namerecord=[\"trans_k\",\"trans_k33\"],\n",
    "    interpolation_methodrecord=[\"linear\",\"linear\"],\n",
    ")\n",
    "\n",
    "npf.write()\n",
    "tvk.write()\n",
    "tvk.ts.write()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af41cba2-bff3-406f-b1a2-dbf87f59621d",
   "metadata": {},
   "source": [
    "# Storage\r\n",
    "The STO package functions in a similar manner to the NPF package but also includes control for stress period type as either transient or steady-state. Instead of conductance behaviour the STO package controls confined versus unconfined behaviour through the 'iconvert' array. To configure specific stress periods for steady-state or transient you need to provide dictionaries keyed to a stress period number with Boolean entries. Note only a change in stress period type needs to be included. The last designated period type remains in effect for all subsequent stress periods unless it is changed\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1eebbdda-5242-4240-b35d-2ea2d0f09835",
   "metadata": {},
   "outputs": [],
   "source": [
    "# building an array for sy\n",
    "sy_layer_prop = [0.01,0.02,0.001]\n",
    "sy_array = np.ones_like(mg.botm)\n",
    "for i,j in enumerate(sy_layer_prop):\n",
    "    sy_array[i]=j\n",
    "\n",
    "# building an ss array\n",
    "ss_array = np.ones_like(mg.botm)*1.0E-5\n",
    "ss_array[1] =  1.0E-6\n",
    "ss_array[2] =  1.0E-7\n",
    "ic_array = np.ones_like(mg.botm)\n",
    "ic_array[2] = 0\n",
    "\n",
    "sto = flopy.mf6.ModflowGwfsto(\n",
    "    gwf,\n",
    "    pname=\"sto\",\n",
    "    save_flows=True,\n",
    "    iconvert=ic_array,\n",
    "    ss=ss_array,\n",
    "    sy=sy_array,\n",
    "    steady_state={0: True},\n",
    "    transient={1: True},\n",
    ")\n",
    "\n",
    "sto.write()\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,f'{model_name}.sto'))]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48f0bb2a-53db-440e-aea0-8f2e3d72c7f2",
   "metadata": {},
   "source": [
    "# Transient storage\n",
    "In keeping with the property change for hydraulic conductivity we can also configure a change to storage with time. This is basically a repeat of what we did previously but specific to the storage package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b18f467-0778-4b5e-8b8e-3c9b1b506192",
   "metadata": {},
   "outputs": [],
   "source": [
    "pdata=[]\n",
    "for lay in range(2):\n",
    "    for node in my_node_idx:\n",
    "        pdata.append(((lay,node),'sy','trans_sy')) # 'sy' is the property I want to change, 'trans_sy' will be the time series that I want to use\n",
    "        pdata.append(((lay,node),'ss','trans_ss')) # 'ss' is the property I want to change, 'trans_ss' will be the time series that I want to use\n",
    "\n",
    "tvs_pdata = {}\n",
    "for key in range(tvk_period,tvk_period+13):\n",
    "    tvs_pdata[key] = pdata\n",
    "tvs = flopy.mf6.ModflowUtltvs(sto, # passing in the sto package\n",
    "                              perioddata=tvs_pdata,\n",
    "                              filename=\"{}.tvs\".format(model_name))\n",
    "\n",
    "# make the tsdata\n",
    "ts_data = [(1.0,0.3,1.0E-04),(stime,0.3,1.0E-04),(sptime+1,0.05,1.0E-06)] # note you must have this start at t=1.0 and end at the end of the model run\n",
    "\n",
    "\n",
    "# initialize first time series\n",
    "tvs.ts.initialize(\n",
    "    filename=\"tvs.ts\",\n",
    "    timeseries=ts_data,\n",
    "    time_series_namerecord=['trans_sy','trans_ss'],\n",
    "    interpolation_methodrecord=[\"linear\",\"linear\"],\n",
    ")\n",
    "\n",
    "sto.write()\n",
    "tvs.write()\n",
    "tvs.ts.write()\n",
    "\n",
    "# take a look at the files written in the model folder and make sure that they are in agreement with the MF6io.pdf document."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7e8103f-3bd8-4be0-b98b-a89b792ffe46",
   "metadata": {},
   "source": [
    "# Recharge Array \r\n",
    "The MF6io.pdf document shows two recharge package types. We demonstrated the \"list type\" in the previous workshop and will now demonstrate the \"array type\". It is useful to know that array type recharge cannot be used with DISU grids. To start off with we are using a mean rainfall value in mm/yr converted to m/d then setting a maximum (3%) and minimum (0.5%) percentage of rainfall as recharge changing over the simulation length to represent a long-term drying trend. Note time array series are different in some aspects to the time series we've used previously. These differences will be highlighted in the comments. This is just one example of what you can do with a combination of multiplier and time series arrays\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "163cc3cf-ea7e-4c82-88f7-561db1dc765c",
   "metadata": {},
   "outputs": [],
   "source": [
    "rain_rate = 1220/365/1000 #m/d\n",
    "\n",
    "aux_period = {} # start a dictionary for an auxiliary mutiplier\n",
    "rch_mult_array=np.ones(np.shape(mg.top)) # create an appropriate array filled with ones\n",
    "aux_period[0]=[rch_mult_array] # set the initial mutiplier to be all 1.0\n",
    "rch_mult_array[my_node_idx] = 3.0 # change the nodes in the zone to have a multiplier of 3.0\n",
    "aux_period[tvk_period] = [rch_mult_array] # set the enhanced recharge to be active from the same period as transient properties\n",
    " \n",
    "rch = flopy.mf6.ModflowGwfrcha( # note the 'a'\n",
    "    gwf,\n",
    "    filename=\"{}.rch\".format(model_name),\n",
    "    pname=\"rch\", # specifying a package name here too\n",
    "    fixed_cell=True,\n",
    "    save_flows=True,\n",
    "    recharge=\"TIMEARRAYSERIES rch_trans\", # Note the use of a keyword followed by the timeseries name record\n",
    "    auxiliary='my_mult_array', # the nmae of our only auxiliary\n",
    "    auxmultname='my_mult_array', # informing the package that our auxiliary is a multiplier\n",
    "    aux=aux_period # passing in the dictionary for our\n",
    ")\n",
    "\n",
    "# creating a time series dictionary keyed to model time and value\n",
    "tas = {0.0: 0.03*rain_rate, # note that instead of a constant value across the whole model this could also be an array with different values\n",
    "       sptime+1: 0.005*rain_rate} # note that instead of a constant value across the whole model this could also be an array with different values\n",
    "\n",
    "rch.tas.initialize(\n",
    "        filename=f\"{model_name}.rch.tas\",\n",
    "        tas_array=tas,\n",
    "        time_series_namerecord=\"rch_trans\", # this must match what you used previously\n",
    "        interpolation_methodrecord=\"LINEAR\",\n",
    ")\n",
    "rch.write()\n",
    "rch.tas.write()\n",
    "# note you can't use observation files with array style recharge and or ET"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bff839fc-8ea1-4579-abdc-e717f4b58288",
   "metadata": {},
   "source": [
    "# ET Array\r\n",
    "Effectively the same as the recharge array package but requires a few extra arrays as input\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd2a9304-164b-4ca5-86cb-85e54faa5a3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "pet = 1600 # mm/yr\n",
    "max_rate = pet/365/1000 # m/d\n",
    "et_rate_array = np.ones_like(mg.top)*max_rate # perhaps a tif file of actual ET could have been used here. \n",
    "\n",
    "ext_depth = 2.5 #meters\n",
    "et_depth_array = np.ones_like(mg.top)*ext_depth # again this could be zoned to represent vegetation maps.\n",
    "\n",
    "evt = flopy.mf6.ModflowGwfevta(\n",
    "    gwf,\n",
    "    readasarrays=True,\n",
    "    fixed_cell=False,\n",
    "    surface = mg.top, # setting the model top as the ET surface\n",
    "    rate = et_rate_array, # passing in our maximum rate array, could have been a timeseries witha  seasonal signal\n",
    "    depth = et_depth_array, # passing our extinction depth array\n",
    "    filename=\"{}.evt\".format(model_name),\n",
    "    pname=\"evt\")\n",
    "evt.write()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "372e7bf4-0516-4b46-b03f-3447228b2ac3",
   "metadata": {},
   "source": [
    "# Initial hydraulic heads\r\n",
    "All numerical models require an estimate of initial hydraulic heads to begin with. A standard approach is to assume top of model to begin with if water tables are not too deep. Using a raster interpolated from measured heads is also easily accomplished using the raster to grid sampling we demonstrated for the surface elevation. Generally, once a steady-state stress period solution is available this commonly be substituted for whatever was used originally as the initial conditions. The code block below demonstrates how to first configure the initial conditions package to leverage the top of model information and then once a heads file is available, use the solution for the steady-state stress period from that file\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "682847df-fefe-471e-b47f-1467c0d0a4ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting up IC array\n",
    "ihd_array=np.ones_like(mg.botm)\n",
    "ihd_array[:] = mg.top-5 # setting the heads in each layer to be 5.0 m below the top of the model, each column of cells feature the same head.\n",
    "\n",
    "# here we are checking the folder that the script is running in for a txt file with initial heads\n",
    "\n",
    "if os.path.exists(\"iheads_array.txt\"):\n",
    "    last_ssheads = np.loadtxt(\"iheads_array.txt\")\n",
    "    print('using previous simualtion SS period heads')\n",
    "    ihd_array[:]=last_ssheads\n",
    "else:\n",
    "    pass\n",
    "\n",
    "ic = flopy.mf6.ModflowGwfic(\n",
    "    gwf, pname=\"ic\", strt=ihd_array, filename=\"{}.ic\".format(model_name)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d09efef-ad06-4396-9aaf-e35d9210363a",
   "metadata": {},
   "source": [
    "# How to get initial heads from a model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc2996d2-a7f3-476e-9eae-295ae3727102",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This will be covered again in a later session on post processing but is included here anyway.\n",
    "# once the model has run to completion you :\n",
    "# 1. create a path to the headsfile\n",
    "# 2. create a hds object from the path\n",
    "# 3. get the data from the stead-state stress period (0,0)\n",
    "# 4. save it as an array with numpy (note we only save the top layer here)\n",
    "\"\"\"\n",
    "headfile = os.path.join(mod_f,\"{}.hds\".format(model_name))\n",
    "hds = flopy.utils.binaryfile.HeadFile(headfile)\n",
    "h = hds.get_data((0,0))\n",
    "np.savetxt(\"iheads_array.txt\",h[0])\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b1f147a-6d98-4ee2-98c4-47e840301da3",
   "metadata": {},
   "source": [
    "# Output control\r\n",
    "The OC package is required to specify what to save and when during your model run. This can be tricky under certain conditions. For the most part you are either going to be saving heads as \"HEAD\" or budgets as \"BUDGET\". In general, you can always save at the start and end of stress periods using keywords \"FIRST\" and \"LAST\". Time step level output is controlled with keywords and/or keyword, integer combinations as will be demonstrated below. Just like most of the other package we created so far, we must build dictionaries with stress period information informing the OC package what to save in each period. Note if you use ATS timing then time step level output may not be possible. Please read the MF6io.pdf document for more information\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b55988d-39b2-4d75-a8b6-10fc45565676",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = list(range(1,numper)) # this range represents the monthly stress periods before a single 100 year stress period\n",
    "hs_keys = [0,*_] # SSkey = 0, then all monthly SP\n",
    "\n",
    "# If I'm not interested in the intercell fluxes then all I need are the heads\n",
    "h_rec = {key:[(\"HEAD\",\"LAST\")] for key in hs_keys} # creating a dictionary to say I want to save the heads at the end of each stress period\n",
    "h_rec[numper] = [(\"HEAD\",\"FREQUENCY\",12)] # the last period which is a 100 year recovery I am using the frequency option for time step level output every 12 time steps\n",
    "\n",
    "# However, if I also want to run something like Zonebudget then I do need the budget so use a combined head plus budget \n",
    "zbud_rec = {key:[(\"BUDGET\",\"LAST\"),(\"HEAD\",\"LAST\")] for key in hs_keys} # note that the entry for each stress period is a list of tuples. \n",
    "zbud_rec[numper] = [(\"BUDGET\",\"FREQUENCY\",12),(\"HEAD\",\"FREQUENCY\",12)] # the numebr of elements in te tuple can vary depending on the options you choose.\n",
    "\n",
    "# A seperate dictionary is used for printing to the list file\n",
    "b_rec = {key:[(\"BUDGET\",\"LAST\")] for key in hs_keys} # this is saying I want a budget summary at the end of each stress period in the listing file.\n",
    "b_rec[numper] = [(\"BUDGET\",\"FREQUENCY\",12)] # this is saying that I also want a budget summary every 12 time steps in the last stress period.  \n",
    "\n",
    "\n",
    "oc = flopy.mf6.ModflowGwfoc(\n",
    "    gwf,\n",
    "    pname=\"oc\",\n",
    "    budget_filerecord=\"{}.cbb\".format(model_name),\n",
    "    head_filerecord=\"{}.hds\".format(model_name),\n",
    "    headprintrecord=[(\"COLUMNS\", 10, \"WIDTH\", 15, \"DIGITS\", 6, \"GENERAL\")], # specify the settings for the a printing record, not actually used here beacuse we are not printing heads to the list file\n",
    "    saverecord=zbud_rec, # here I am using the dictionary that will also save a budget file\n",
    "    printrecord=b_rec,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0844cbd-a4ab-4399-b385-53f336264725",
   "metadata": {},
   "source": [
    "# Model observations\r\n",
    "We demonstrated how to create stress package observation files but there is also a model level observation package. These are useful for obtaining heads at specific model cells for hydrograph plotting. You can also get drawdown as temporal differences or net lateral flux from a cell. The model level observations are either heads, flux or drawdown from specific cells\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54e47193-e211-44fa-9d9b-943445656981",
   "metadata": {},
   "outputs": [],
   "source": [
    "# using intersection with some geopandas to get my observation model cells.\n",
    "ix = GridIntersect(mg, method=\"vertex\")\n",
    "my_shape = os.path.join(gis_f,'L1_monitoring.shp') # a point shapefile with well locations\n",
    "gdf = gpd.read_file(my_shape)\n",
    "my_poly = gpd.read_file(my_shape).geometry # note with a point shapefile we don't pick the zero index like we did previously\n",
    "my_cells = [ix.intersect(x) for x in my_poly]\n",
    "my_node_idx = [x[0][0] for x in my_cells] # but here we have an extra zero index to get just the node indices as a list\n",
    "\n",
    "# creating the observation data list from the modelgrid intersection information\n",
    "names = gdf['Name'].tolist()\n",
    "nodes = my_node_idx\n",
    "assert len(nodes)==len(names) # making sure that my nodes and names are equal\n",
    "obs_list = [(names[i],\"HEAD\",1,nodes[i]) for i in range(len(names))] # for some wierd reason the layer number needs to not be zero based.\n",
    "\n",
    "# This package is built the same as most others.\n",
    "obs_dict = {} \n",
    "obs_dict['head_obs.csv'] = obs_list\n",
    "obs_package = flopy.mf6.ModflowUtlobs(\n",
    "     gwf,\n",
    "     pname=\"head_obs\",\n",
    "     filename=\"{}.obs\".format(model_name),\n",
    "     digits=10,\n",
    "     print_input=True,\n",
    "     continuous=obs_dict,\n",
    " )\n",
    "obs_package.write()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a21d2853-ea18-4e87-8fdd-151f78ff8a1e",
   "metadata": {},
   "source": [
    "# Basics covered\n",
    "That should be all you really need to know how to do in order to get your models up and running. \n",
    "Our plan for the next sessions are as follows:\n",
    "1. Session 5 = Advanced stress packages plus basic post processing\n",
    "2. Session 6 = Solute and variable density transport\n",
    "3. Session 7 = ZoneBudget and MODPATH\n",
    "4. Session 8 = Putting all together\n",
    "\n",
    "See you next week."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "159a1f7c-d004-41d4-bf7b-5848808f9492",
   "metadata": {},
   "outputs": [],
   "source": [
    "shutil.rmtree(ws4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5561ceb1-9ded-4130-8d65-157f5d6bd741",
   "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
}