examples/Notebooks/flopy3_mt3d-usgs_example_with_sft_lkt_uzt.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demonstrate use of SFR/LAK/UZF with SFT/LKT/UZT\n",
    "A more comprehensive demonstration of setting up an MT3D-USGS model that uses all of the new packages included in the first release of MT3D-USGS.  Also includes RCT.\n",
    "\n",
    "#### Problem Description:\n",
    "* 300 row x 300 col x 3 layer x 2 stress period model\n",
    "* Flow model uses SFR, LAK, and UZF with connections between all three\n",
    "* Transport model simulates streamflow transport (SFT), with connction to a single lake (LKT)\n",
    "* Transport model simulates overland runoff and spring discharge (UZT) to surface water network\n",
    "\n",
    "Start by importing some libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.8.6 | packaged by conda-forge | (default, Oct  7 2020, 18:42:56) \n",
      "[Clang 10.0.1 ]\n",
      "numpy version: 1.18.5\n",
      "matplotlib version: 3.2.2\n",
      "flopy version: 3.3.3\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "import platform\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# run installed version of flopy or add local path\n",
    "try:\n",
    "    import flopy\n",
    "except:\n",
    "    fpth = os.path.abspath(os.path.join('..', '..'))\n",
    "    sys.path.append(fpth)\n",
    "    import flopy\n",
    "\n",
    "print(sys.version)\n",
    "print('numpy version: {}'.format(np.__version__))\n",
    "print('matplotlib version: {}'.format(mpl.__version__))\n",
    "print('flopy version: {}'.format(flopy.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a MODFLOW model and store it, in this case in the variable 'mf'. \n",
    "The modelname will be the name given to all MODFLOW files.\n",
    "The exe_name should be the name of the MODFLOW executable. \n",
    "In this case, we want to use version: 'mfnwt' for MODFLOW-NWT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "modelpth = os.path.join('.', 'temp', 'no3')\n",
    "modelname = 'no3'\n",
    "mfexe = 'mfnwt'\n",
    "mtexe = 'mt3dusgs'\n",
    "\n",
    "if platform.system() == 'Windows':\n",
    "    mfexe += '.exe'\n",
    "    mtexe += '.exe'\n",
    "\n",
    "# Make sure modelpth directory exists\n",
    "if not os.path.exists(modelpth):\n",
    "    os.makedirs(modelpth)\n",
    "\n",
    "# Instantiate MODFLOW object in flopy\n",
    "mf = flopy.modflow.Modflow(modelname=modelname, exe_name=mfexe, model_ws=modelpth, version='mfnwt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up model discretization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lx = 90000.0 \n",
    "Ly = 90000.0\n",
    "nrow = 300\n",
    "ncol = 300\n",
    "nlay = 3\n",
    "\n",
    "delr = Lx / ncol\n",
    "delc = Ly / nrow\n",
    "\n",
    "xmax = ncol * delr\n",
    "ymax = nrow * delc\n",
    "\n",
    "X, Y = np.meshgrid(np.linspace(delr / 2, xmax - delr / 2, ncol),    \n",
    "                   np.linspace(ymax - delc / 2, 0 + delc / 2, nrow))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate output control (oc) package for MODFLOW-NWT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "oc = flopy.modflow.ModflowOc(mf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate solver package for MODFLOW-NWT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Newton-Raphson Solver: Create a flopy nwt package object\n",
    "\n",
    "headtol = 1.0E-4 \n",
    "fluxtol = 5      \n",
    "maxiterout = 5000\n",
    "thickfact = 1E-06\n",
    "linmeth = 2      \n",
    "iprnwt = 1       \n",
    "ibotav = 1       \n",
    "\n",
    "nwt = flopy.modflow.ModflowNwt(mf, headtol=headtol, fluxtol=fluxtol, maxiterout=maxiterout,\n",
    "                               thickfact=thickfact, linmeth=linmeth, iprnwt=iprnwt, ibotav=ibotav,\n",
    "                               options='SIMPLE')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate discretization (DIS) package for MODFLOW-NWT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "elv_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','dis_arrays','grnd_elv.txt')\n",
    "\n",
    "# Top of Layer 1 elevation determined using GW Vistas and stored locally\n",
    "grndElv = np.loadtxt(elv_pth)\n",
    "\n",
    "# Bottom of layer 1 elevation also determined from use of GUI and stored locally\n",
    "bt1_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','dis_arrays','bot1.txt')\n",
    "bot1Elv = np.loadtxt(bt1_pth)\n",
    "\n",
    "bot2Elv = np.ones(bot1Elv.shape) * 100\n",
    "bot3Elv = np.zeros(bot2Elv.shape)\n",
    "\n",
    "botm = [bot1Elv, bot2Elv, bot3Elv]\n",
    "botm = np.array(botm)\n",
    "Steady = [False, False]\n",
    "nstp = [1, 1]\n",
    "tsmult = [1., 1.]\n",
    "\n",
    "# Stress periods \n",
    "perlen = [9131.25, 9131.25]\n",
    "\n",
    "# Create the discretization object\n",
    "# itmuni = 4 (days); lenuni = 1 (feet)\n",
    "dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, nper=2, delr=delr, delc=delc,     \n",
    "                               top=grndElv, botm=botm, laycbd=0, itmuni=4, lenuni=1,\n",
    "                               steady=Steady, nstp=nstp, tsmult=tsmult, perlen=perlen) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate upstream weighting (UPW) flow package for MODFLOW-NWT\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# UPW must be instantiated after DIS.  Otherwise, during the mf.write_input() procedures,   \n",
    "# flopy will crash.\n",
    "\n",
    "# First line of UPW input is: IUPWCB HDRY NPUPW IPHDRY \n",
    "hdry = -1.00e+30\n",
    "iphdry = 0\n",
    "\n",
    "# Next variables are: LAYTYP, LAYAVG, CHANI, LAYVKA, LAYWET\n",
    "laytyp = [1, 3, 3]  # >0: convertible\n",
    "layavg = 0          #  0: harmonic mean\n",
    "chani = 1.0         # >0: CHANI is the horizontal anisotropy for the entire layer\n",
    "layvka = 0          # =0: indicates VKA is vertical hydraulic conductivity\n",
    "laywet = 0          # Always set equal to zero in UPW package\n",
    "hk = 20\n",
    "# hani = 1          # Not needed because CHANI > 1\n",
    "vka = 0.5           # Is equal to vert. K b/c LAYVKA = 0\n",
    "ss = 0.00001\n",
    "sy = 0.20\n",
    "\n",
    "upw = flopy.modflow.ModflowUpw(mf, laytyp=laytyp, layavg=layavg, chani=chani, layvka=layvka,\n",
    "                               laywet=laywet, ipakcb=53, hdry=hdry, iphdry=iphdry, hk=hk, \n",
    "                               vka=vka, ss=ss, sy=sy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate basic (BAS or BA6) package for MODFLOW-NWT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "ibnd1_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','bas_arrays','ibnd_lay1.txt')\n",
    "ibnd1 = np.loadtxt(ibnd1_pth)\n",
    "ibnd2 = np.ones(ibnd1.shape)\n",
    "ibnd3 = np.ones(ibnd2.shape)\n",
    "\n",
    "ibnd = [ibnd1, ibnd2, ibnd3]\n",
    "ibnd = np.array(ibnd)\n",
    "\n",
    "StHd1_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','bas_arrays','strthd1.txt')\n",
    "StHd1 = np.loadtxt(StHd1_pth)\n",
    "\n",
    "StHd2_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','bas_arrays','strthd2.txt')\n",
    "StHd2 = np.loadtxt(StHd2_pth)\n",
    "\n",
    "StHd3_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','bas_arrays','strthd3.txt')\n",
    "StHd3 = np.loadtxt(StHd3_pth)\n",
    "\n",
    "strtElev = [StHd1, StHd2, StHd3]\n",
    "strtElev = np.array(strtElev)\n",
    "\n",
    "hdry = 999.0\n",
    "\n",
    "bas = flopy.modflow.ModflowBas(mf, ibound=ibnd, hnoflo=hdry, strt=strtElev)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate general head boundary (GHB) package for MODFLOW-NWT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GHB boundaries are located along the top (north) and bottom (south) \n",
    "# edges of the domain, all 3 layers.\n",
    "\n",
    "elev_stpt_row1   = 308.82281\n",
    "elev_stpt_row300 = 239.13811\n",
    "elev_slp         = (308.82281 - 298.83649) / (ncol - 1)\n",
    "\n",
    "sp = []\n",
    "for k in [0,1,2]:                     # These indices need to be adjusted for 0-based moronicism\n",
    "    for i in [0,299]:                 # These indices need to be adjusted for 0-based silliness\n",
    "        for j in np.arange(0,300,1):  # These indices need to be adjusted for 0-based foolishness\n",
    "            # Skipping cells not satisfying the conditions below\n",
    "            if ((i == 1 and (j < 27 or j > 31)) or (i==299 and (j < 26 or j > 31))):\n",
    "                if (i % 2 == 0):\n",
    "                    sp.append([k, i, j, elev_stpt_row1 - (elev_slp * (j-1)), 11.3636])\n",
    "                else:\n",
    "                    sp.append([k, i, j, elev_stpt_row300 - (elev_slp * (j-1)), 11.3636])\n",
    "\n",
    "\n",
    "for k in [0,1,2]:\n",
    "    for j in np.arange(26,32,1):\n",
    "        sp.append([k, 299, j, 238.20, 3409.0801])\n",
    "\n",
    "ghb = flopy.modflow.ModflowGhb(mf, stress_period_data = sp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate streamflow routing (SFR2) package for MODFLOW-NWT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read pre-prepared reach data into numpy recarrays using numpy.genfromtxt()\n",
    "# Remember that the cell indices stored in the pre-prepared NO3_ReachInput.csv file are based on 0-based indexing.\n",
    "# Flopy will convert to 1-based when it writes the files\n",
    "\n",
    "rpth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','sfr_data', 'no3_reachinput.csv')\n",
    "reach_data = np.genfromtxt(rpth, delimiter=',', names=True)\n",
    "reach_data\n",
    "\n",
    "# Read pre-prepared segment data into numpy recarrays using numpy.genfromtxt()\n",
    "\n",
    "spth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','sfr_data', 'no3_segmentdata.csv')\n",
    "ss_segment_data = np.genfromtxt(spth, delimiter=',', names=True)\n",
    "segment_data = {0: ss_segment_data,\n",
    "                1: ss_segment_data}\n",
    "segment_data[0][0:1]['width1']\n",
    "\n",
    "nstrm = len(reach_data)\n",
    "nss = len(segment_data[0])\n",
    "nsfrpar = 0\n",
    "const = 128390.4  # constant for manning's equation, units of cfs\n",
    "dleak = 0.0001\n",
    "ipakcb = 53       # flag for writing SFR output to cell-by-cell budget (on unit 53)\n",
    "istcb2 = 37       # flag for writing SFR output to text file\n",
    "isfropt = 1\n",
    "dataset_5 = {0: [nss, 0, 0], 1: [-1, 0, 0]} # dataset 5 (see online guide) (ITMP, IRDFLG, IPTFLG)\n",
    "\n",
    "# Input arguments generally follow the variable names defined in the Online Guide to MODFLOW\n",
    "sfr = flopy.modflow.ModflowSfr2(mf, nstrm=nstrm, nss=nss, const=const, dleak=dleak, ipakcb=ipakcb,\n",
    "                                istcb2=istcb2, isfropt=isfropt, reachinput=True,\n",
    "                                reach_data=reach_data,\n",
    "                                segment_data=segment_data,\n",
    "                                dataset_5=dataset_5, unit_number=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate Lake (LAK) package for MODFLOW-NWT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read pre-prepared lake arrays\n",
    "LakArr_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','lak_arrays','lakarr1.txt')\n",
    "LakArr_lyr1 = np.loadtxt(LakArr_pth)\n",
    "LakArr_lyr2 = np.zeros(LakArr_lyr1.shape)\n",
    "LakArr_lyr3 = np.zeros(LakArr_lyr2.shape)\n",
    "\n",
    "LakArr = [LakArr_lyr1, LakArr_lyr2, LakArr_lyr3]\n",
    "LakArr = np.array(LakArr)\n",
    "\n",
    "nlakes  = int(np.max(LakArr))\n",
    "ipakcb  = ipakcb      # From above\n",
    "theta   = -1.         # Implicit\n",
    "nssitr  = 10          # Maximum number of iterations for Newton’s method\n",
    "sscncr  = 1.000e-03   # Convergence criterion for equilibrium lake stage solution\n",
    "surfdep = 2.000e+00   # Height of small topological variations in lake-bottom\n",
    "stages  = 268.00      # Initial stage of each lake at the beginning of the run\n",
    "\n",
    "# ITMP  > 0, read lake definition data \n",
    "# ITMP1 ≥ 0, read new recharge, evaporation, runoff, and withdrawal data for each lake\n",
    "# LWRT  > 0, suppresses printout from the lake package\n",
    "\n",
    "bdlknc_lyr1 = LakArr_lyr1.copy()\n",
    "bdlknc_lyr2 = LakArr_lyr1.copy()\n",
    "bdlknc_lyr3 = np.zeros((LakArr_lyr1.shape))\n",
    "\n",
    "# Need to expand bdlknc_lyr1 non-zero values by 1 in either direction\n",
    "# (left/right and up/down)\n",
    "for i in np.arange(0, LakArr_lyr1.shape[0]):\n",
    "    for j in np.arange(0, LakArr_lyr1.shape[1]):\n",
    "        im1 = i - 1\n",
    "        ip1 = i + 1\n",
    "        jm1 = j - 1\n",
    "        jp1 = j + 1\n",
    "        \n",
    "        if(im1 >= 0):\n",
    "            if(LakArr_lyr1[i, j] == 1 and LakArr_lyr1[im1, j] == 0):\n",
    "                bdlknc_lyr1[im1, j] = 1\n",
    "        \n",
    "        if(ip1 < LakArr_lyr1.shape[0]):\n",
    "            if(LakArr_lyr1[i, j] == 1 and LakArr_lyr1[ip1, j] == 0):\n",
    "                bdlknc_lyr1[ip1, j] = 1\n",
    "        \n",
    "        if(jm1 >= 0):\n",
    "            if(LakArr_lyr1[i, j] == 1 and LakArr_lyr1[i, jm1] == 0):\n",
    "                bdlknc_lyr1[i, jm1] = 1\n",
    "        \n",
    "        if(jp1 < LakArr_lyr1.shape[1]):\n",
    "            if(LakArr_lyr1[i, j] == 1 and LakArr_lyr1[i, jp1] == 0):\n",
    "                bdlknc_lyr1[i, jp1] = 1\n",
    "\n",
    "\n",
    "bdlknc = [bdlknc_lyr1, bdlknc_lyr2, bdlknc_lyr3]\n",
    "bdlknc = np.array(bdlknc)\n",
    "\n",
    "flux_data = {0: [[0.0073, 0.0073, 0.0, 0.0]],\n",
    "             1: [[0.0073, 0.0073, 0.0, 0.0]]}\n",
    "\n",
    "lak = flopy.modflow.ModflowLak(mf, nlakes=nlakes, ipakcb=ipakcb, theta=theta,\n",
    "                                nssitr=nssitr, sscncr=sscncr, surfdep=surfdep,\n",
    "                                stages=stages, lakarr=LakArr, bdlknc=bdlknc,\n",
    "                                flux_data=flux_data, unit_number=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate gage package for use with MODFLOW-NWT package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "gages = [[1, 225, 90, 3],\n",
    "         [2,  68, 91, 3],\n",
    "         [3,  33, 92, 3],\n",
    "         [4, 165, 93, 3],\n",
    "         [5, 123, 94, 3],\n",
    "         [6,  77, 95, 3],\n",
    "         [7, 173, 96, 3],\n",
    "         [8, 328, 97, 3],\n",
    "         [9, 115, 98, 3],\n",
    "         [-1,   -101, 1]]\n",
    "\n",
    "# gages = [[1,38,61,1],[2,67,62,1], [3,176,63,1], [4,152,64,1], [5,186,65,1], [6,31,66,1]]\n",
    "files = ['no3.gag','seg1_gag.out','seg2_gag.out','seg3_gag.out','seg4_gag.out','seg5_gag.out','seg6_gag.out',\n",
    "         'seg7_gag.out','seg8_gag.out','seg9_gag.out','lak1_gag.out']\n",
    "\n",
    "numgage = len(gages)\n",
    "gage = flopy.modflow.ModflowGage(mf, numgage=numgage, gage_data=gages, filenames = files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate Unsaturated-Zone Flow (UZF) package for MODFLOW-NWT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "nuztop  =  2\n",
    "iuzfopt =  2\n",
    "irunflg =  1\n",
    "ietflg  =  0\n",
    "iuzfcb  = 52\n",
    "iuzfcb2 =  0\n",
    "ntrail2 = 20\n",
    "nsets2  = 20\n",
    "nuzgag  =  2\n",
    "surfdep =  2.\n",
    "\n",
    "eps     =  3.\n",
    "thts    =  0.30\n",
    "thti    =  0.13079\n",
    "\n",
    "fname_uzbnd = os.path.join('..','data','mt3d_example_sft_lkt_uzt','uzf_arrays','iuzbnd.txt')\n",
    "fname_runbnd = os.path.join('..','data','mt3d_example_sft_lkt_uzt','uzf_arrays','irunbnd.txt')\n",
    "\n",
    "iuzfbnd = np.loadtxt(fname_uzbnd)\n",
    "irunbnd = np.loadtxt(fname_runbnd)\n",
    "\n",
    "uzgag = [[106, 160,  121, 3],\n",
    "         [  1,   1, -122, 1]]\n",
    "\n",
    "finf = {0: 1.8250e-03, 1: 1.8250e-03}\n",
    "\n",
    "uzf = flopy.modflow.ModflowUzf1(mf, nuztop=nuztop, iuzfopt=iuzfopt, irunflg=irunflg, ietflg=ietflg, ipakcb=iuzfcb,\n",
    "                               iuzfcb2=iuzfcb2, ntrail2=ntrail2, nsets=nsets2, surfdep=surfdep, uzgag=uzgag,\n",
    "                               iuzfbnd=1, irunbnd=0, vks=1.0E-6, eps=3.5, thts=0.35,\n",
    "                               thtr=0.15, thti=0.20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate Drain (DRN) package for MODFLOW-NWT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname_drnElv = os.path.join('..','data','mt3d_example_sft_lkt_uzt','drn_arrays', 'elv.txt')\n",
    "fname_drnCond = os.path.join('..','data','mt3d_example_sft_lkt_uzt','drn_arrays','cond.txt')\n",
    "\n",
    "drnElv = np.loadtxt(fname_drnElv)\n",
    "drnCond = np.loadtxt(fname_drnCond)\n",
    "\n",
    "drnElv_lst = pd.DataFrame({\n",
    "                 'lay':  1,\n",
    "                 'row':  np.nonzero(drnElv)[0] + 1,\n",
    "                 'col':  np.nonzero(drnElv)[1] + 1,\n",
    "                 'elv':  drnElv[np.nonzero(drnElv)],\n",
    "                 'cond': drnCond[np.nonzero(drnCond)] }, \n",
    "             columns=['lay', 'row', 'col', 'elv', 'cond'])\n",
    "\n",
    "# Convert the DataFrame into a list of lists for the drn constructor\n",
    "stress_period_data = drnElv_lst.values.tolist()\n",
    "\n",
    "# Create a dictionary, 1 entry for each of the two stress periods.\n",
    "stress_period_data = {0: stress_period_data,\n",
    "                      1: stress_period_data}\n",
    "\n",
    "drn = flopy.modflow.ModflowDrn(mf, ipakcb=ipakcb, stress_period_data=stress_period_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate linkage with mass transport routing (LMT) package for MODFLOW-NWT (generates linker file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "lmt = flopy.modflow.ModflowLmt(mf, output_file_name='NO3.ftl', output_file_header='extended',\n",
    "                               output_file_format='formatted', package_flows = ['all'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now work on MT3D-USGS file creation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start by setting up MT3D-USGS class and pass in MODFLOW-NWT object for setting up a number of the BTN arrays\n",
    "\n",
    "mt = flopy.mt3d.Mt3dms(modflowmodel=mf, modelname=modelname, model_ws=modelpth,\n",
    "                       version='mt3d-usgs', namefile_ext='mtnam', exe_name=mtexe,\n",
    "                       ftlfilename='no3.ftl', ftlfree=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate basic transport (BTN) package for MT3D-USGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncomp   =   1\n",
    "lunit   =   'FT'\n",
    "sconc   =   0.0\n",
    "prsity  =   0.3\n",
    "cinact  =  -1.0\n",
    "thkmin  =   0.000001\n",
    "nprs    =  -2\n",
    "nprobs  =  10\n",
    "nprmas  =  10\n",
    "dt0     =   0.1\n",
    "nstp    =   1\n",
    "mxstrn  = 500\n",
    "ttsmult =   1.2\n",
    "ttsmax  = 100\n",
    "\n",
    "# These observations need to be entered with 0-based indexing\n",
    "obs = [[0, 104, 158],\n",
    "       [1, 104, 158],\n",
    "       [2, 104, 158]]\n",
    "\n",
    "btn = flopy.mt3d.Mt3dBtn(mt, MFStyleArr=True, DRYCell=True, lunit=lunit,\n",
    "                         sconc=sconc, ncomp=ncomp, prsity=prsity, cinact=cinact,\n",
    "                         obs=obs, thkmin=thkmin, nprs=nprs, nprobs=nprobs, \n",
    "                         chkmas=True,nprmas=nprmas, dt0=dt0, nstp=nstp, \n",
    "                         mxstrn=mxstrn, ttsmult=ttsmult, ttsmax=ttsmax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate advection (ADV) package for MT3D-USGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "mixelm =    0\n",
    "percel =    1.0000\n",
    "mxpart = 5000\n",
    "nadvfd =    1      # (1 = Upstream weighting)\n",
    "\n",
    "adv = flopy.mt3d.Mt3dAdv(mt, mixelm=mixelm, percel=percel, mxpart=mxpart, \n",
    "                         nadvfd=nadvfd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate generalized conjugate gradient solver (GCG) package for MT3D-USGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "mxiter =  1\n",
    "iter1  = 50\n",
    "isolve =  3\n",
    "ncrs   =  0\n",
    "accl   =  1.000000\n",
    "cclose =  1.00e-06\n",
    "iprgcg =  5\n",
    "\n",
    "gcg = flopy.mt3d.Mt3dGcg(mt, mxiter=mxiter, iter1=iter1, \n",
    "                         isolve=isolve, ncrs=ncrs, accl=accl, \n",
    "                         cclose=cclose, iprgcg=iprgcg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate dispersion (DSP) package for MT3D-USGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "al     = 0.1  # longitudinal dispersivity\n",
    "trpt   = 0.1  # ratio of the horizontal transverse dispersivity to 'AL'\n",
    "trpv   = 0.1  # ratio of the vertical transverse dispersitvity to 'AL'\n",
    "dmcoef = 1.0000e-10\n",
    "\n",
    "dsp = flopy.mt3d.Mt3dDsp(mt, al=al, trpt=trpt, trpv=trpv, dmcoef=dmcoef,\n",
    "                         multiDiff=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate source-sink mixing (SSM) package for MT3D-USGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# no user-specified concentrations associated with boundary conditions\n",
    "\n",
    "mxss = 11199\n",
    "\n",
    "ssm = flopy.mt3d.Mt3dSsm(mt, mxss=mxss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate reaction (RCT) package for MT3D-USGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "isothm    = 0\n",
    "ireact    = 1\n",
    "irctop    = 2\n",
    "igetsc    = 0\n",
    "ireaction = 0\n",
    "\n",
    "rc1       = 6.3258e-04  # first-order reaction rate for the dissolved phase\n",
    "rc2       = 0.0         # Decay on Soil Layer\n",
    "\n",
    "rct = flopy.mt3d.Mt3dRct(mt, isothm=isothm, ireact=ireact, igetsc=igetsc,\n",
    "                         rc1=rc1, rc2=rc2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate streamflow transport (SFT) package for MT3D-USGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsfinit  = len(reach_data)\n",
    "mxsfbc   = len(reach_data)\n",
    "icbcsf   = 0\n",
    "ioutobs  = 92\n",
    "isfsolv  = 1\n",
    "wimp     = 0.5\n",
    "wups     = 1.0\n",
    "cclosesf = 1.0E-6\n",
    "mxitersf = 10\n",
    "crntsf   = 1.0\n",
    "iprtxmd  = 0\n",
    "coldsf   = 0\n",
    "dispsf   = 0\n",
    "obs_sf   = [225,293,326,491,614,691,864,1192,1307]\n",
    "sf_stress_period_data = {0: [0, 0, 0],\n",
    "                         1: [0, 0, 0],\n",
    "                         2: [0, 0, 0]}\n",
    "\n",
    "gage_output = [None, None, 'no3.sftobs']\n",
    "\n",
    "sft = flopy.mt3d.Mt3dSft(mt, nsfinit=nsfinit, mxsfbc=mxsfbc, icbcsf=icbcsf,\n",
    "                         ioutobs=ioutobs, isfsolv=isfsolv, wimp=wimp,\n",
    "                         wups=wups, cclosesf=cclosesf, mxitersf=mxitersf,\n",
    "                         crntsf=crntsf, iprtxmd=iprtxmd, coldsf=coldsf,\n",
    "                         dispsf=dispsf, nobssf=len(obs_sf), obs_sf=obs_sf,\n",
    "                         sf_stress_period_data = sf_stress_period_data,\n",
    "                         filenames=gage_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate unsaturated-zone transport (UZT) package for MT3D-USGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "mxuzcon = np.count_nonzero(irunbnd)\n",
    "icbcuz  = 45\n",
    "iet     = 0\n",
    "wc      = np.ones((nlay, nrow, ncol)) * 0.29\n",
    "sdh     = np.ones((nlay, nrow, ncol))\n",
    "\n",
    "uzt = flopy.mt3d.Mt3dUzt(mt, mxuzcon=mxuzcon, icbcuz=icbcuz, iet=iet,\n",
    "                         iuzfbnd=iuzfbnd, sdh=sdh, cuzinf=1.4158e-03,\n",
    "                         filenames='no3')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate lake transport (LKT) package for MT3D-USGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "nlkinit = 1\n",
    "mxlkbc  = 720\n",
    "icbclk  = 81\n",
    "ietlak  = 1\n",
    "coldlak = 1\n",
    "\n",
    "lkt_flux_data = {0: [[0, 1, 0.01667]],\n",
    "                 1: [[0, 1, 0.02667]]}\n",
    "\n",
    "lkt = flopy.mt3d.Mt3dLkt(mt, nlkinit=nlkinit, mxlkbc=mxlkbc, icbclk=icbclk,\n",
    "                         ietlak=ietlak, coldlak=coldlak, \n",
    "                         lk_stress_period_data=lkt_flux_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Write the MT3D-USGS input files for inspecting and running"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/jdhughes/Documents/Development/flopy_git/flopy_fork/examples/Notebooks\n"
     ]
    }
   ],
   "source": [
    "pth = os.getcwd()\n",
    "print(pth)\n",
    "\n",
    "mf.write_input()\n",
    "mt.write_input()\n",
    "\n",
    "#mf.run_model()\n",
    "#mt.run_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}